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PREFACE 


The  objective  of  this  research  project  was  to  develop  different  geophysical  techniques  for 
environmental  site  characterization.  The  effort  was  concentrated  in  four  areas:  (1)  a  multi¬ 
offset  ground  penetrating  radar  (GPR)  profiling  technique;  (2)  a  nonlinear  seismic  refraction 
traveltime  tomography;  (3)  a  3-D  d.c.  resistivity  tomography;  and  (4)  an  electroseismic 
method.  The  goals  of  this  project  were  accomplished  successfully.  In  addition,  a  large 
amount  of  field  data  were  collected  and  modeled  for  characterizing  actual  environmental 
sites. 

In  the  first  section,  we  demonstrate  that  performing  GPR  surveys  with  multi-offset  ge¬ 
ometry  can  improve  the  signal-to-noise  ratio  of  subsurface  radar  reflections.  Furthermore, 
we  show  that  using  a  common  midpoint  (CMP)  processing  technique  to  model  multi-offset 
radar  data  allows  for  the  determination  of  dielectric  properties,  an  interpretation  of  subsur¬ 
face  water  content,  and  the  conversion  of  radar  traveltime  profiles  to  depth  profiles. 

In  the  second  section,  we  formulate  a  nonlinear  seismic  refraction  traveltime  tomography 
to  precisely  delineate  bedrock  topography  that  fits  not  only  the  traveltime  data  but  also  the 
gradients  of  traveltime  curves.  The  use  of  the  Tikhonov  regularization  method  ensures  that 
the  inverse  problem  is  well-posed  for  models  with  a  large  number  of  grids. 

The  third  section  describes  a  3-D  d.c.  resistivity  tomography  technique  for  characteriz¬ 
ing  conductivity  distribution  in  the  shallow  earth  that  can  map  surface  measurements  to  the 
subsurface.  We  investigate  the  performance  of  different  inversion  algorithms  for  rapidly  mod- 
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eling  a  large  number  of  data  points  on  a  workstation  and  on  an  nCUBE  parallel  computer. 
A  nonlinear  conjugate  gradient  method  performs  consistently  well. 

Finally,  the  fourth  section  of  the  report  is  on  the  theoretical  and  numerical  techniques  for 
modeling  electroseismic  waves  from  a  point  source  by  solving  coupled  equations  of  seismic 
wave  propagation,  fluid  flow,  electrical  streaming  current,  and  the  resulting  electromagnetic 
field. 

The  usefulness  of  these  techniques  for  characterizing  environmental  sites  is  demonstrated 
by  applications  to  field  data  from  various  sites.  The  new  approaches  and  algorithms  intro¬ 
duced  for  GPR,  seismic  refraction  tomography,  and  3-D  resistivity  inversion  have  received 
wide  interest  from  both  the  scientific  and  user  communities.  The  electroseismic  method  is 
still  a  new  and  potentially  powerful  technique  for  subsurface  characterization.  The  theoret¬ 
ical  modeling  that  was  done  is  the  first  of  its  kind  in  this  field. 


viii 


VELOCITY  VARIATIONS  AND  WATER  CONTENT 
ESTIMATED  FROM  MULTI-OFFSET  GROUND 
PENETRATING  RADAR 


Summary 

The  common  midpoint  (CMP)  processing  technique  has  been  shown  to  be  effective  in  im¬ 
proving  the  results  of  ground  penetrating  radar  profiling.  When  radar  data  are  collected 
with  the  CMP  multi-offset  geometry,  stacking  increases  the  signal-to-noise  ratio  of  subsur¬ 
face  radar  reflections  and  results  in  an  improved  subsurface  image.  An  important  aspect  of 
CMP  processing  is  normal  moveout  velocity  analysis.  Our  objectives  are  to  show  the  effect 
of  multiple  velocity  analyses  on  the  stacked  radar  image  and  particularly,  to  demonstrate 
that  this  velocity  information  can  also  be  used  to  determine  subsurface  water  content. 

Most  GPR  surveys  are  very  limited  in  spatial  extent  and  assume  that  within  the  survey 
range,  radar  velocity  structure  in  the  shallow  subsurface  can  be  adequately  approximated 
by  a  single  velocity  function  in  data  processing.  In  this  study  we  show  that  variation  in 
radar  velocity  can  be  quite  significant  and  that  the  stacked  profile  improves  as  the  number 
of  velocity  analysis  locations  is  increased. 

Interval  velocities  can  be  calculated  from  the  normal  moveout  velocities  derived  in  the 
CMP  velocity  analysis.  With  some  reasonable  assumptions  about  subsurface  conditions  nec¬ 
essary  for  radar  propagation,  interval  velocity  can  be  converted  to  an  estimate  of  volumetric 
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water  content.  Therefore,  by  collecting  GPR  data  in  the  multi-offset  CMP  geometry,  not 
only  is  the  radar  profile  improved  but  it  also  allows  for  an  interpretation  of  subsurface  vari¬ 
ation  in  water  content.  We  show  the  application  of  these  techniques  to  multi-offset  GPR 
data  from  the  Chalk  River  test  area  operated  by  Atomic  Energy  of  Canada  Limited. 

1.  Introduction 

Ground  penetrating  radar  (GPR)  is  the  recording  of  high  frequency  electromagnetic  waves 
that  have  reflected  from  subsurface  contrasts  in  dielectric  constant.  In  low-loss  media,  ground 
penetrating  radar  is  capable  of  producing  high-resolution  images  of  the  shallow  subsurface, 
compared  to  other  surface  electric  or  seismic  geophysical  techniques.  Most  commonly,  GPR 
data  is  collected  with  a  constant  source-to-receiver  oflfeet.  In  the  early  development  of  GPR, 
multi-offset  common  midpoint  (CMP)  sounding  was  borrowed  from  reflection  seismology  as 
an  effective  technique  for  determining  glacial  ice  velocities  (Gudmandsen,  1971).  As  the  use 
of  GPR  expanded  to  rock  and  soil  surveys,  multi-offset  radar  sounding  has  continued  to  be 
used  primarily  for  velocity  sounding  at  one  or  just  a  few  points  along  a  survey  hue.  The 
advantages  of  determining  velocity  with  this  method  are  that  it  requires  no  prior  knowledge 
of  the  subsurface,  is  not  intrusive,  uses  the  radar  data  acquisition  system  only  and  can 
determine  a  velocity  profile  anywhere  within  the  survey.  The  disadvantage  is  that  acquiring 
multi-offset  data  with  current  GPR  systems  is  slow  (compared  to  constant  offset  surveys), 
or  impossible,  for  systems  with  fixed  source-receiver  offset. 
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Recent  case  studies  have  shown  that  when  an  entire  GPR  survey  is  acquired  with  the 
CMP  geometry,  multi-trace  reflection  seismic  processing  techniques  can  be  used  to  improve 
the  subsurface  radar  images  (Fisher  et  al,  1992;  Gerlitz  et  al,  1993).  Fisher  et  al.  showed 
the  improvement  of  the  subsurface  image  resulting  from  CMP  stacking  and  migration.  Our 
focus  is  on  utilization  of  the  subsurface  velocity  determined  from  the  CMP  pre-stack  data. 
When  the  data  are  collected  in  the  CMP  configuration,  normal  moveout  velocity  analysis 
is  used  to  derive  a  2-D  radar  stacking  velocity  field.  In  our  study,  using  the  same  data  set 
as  Fisher  et  al.,  we  show  how  the  subsurface  image,  resulting  from  the  stacking  process, 
improves  as  the  number  of  normal  moveout  velocity  analyses,  defining  the  stacking  velocity 
field,  increases.  We  then  demonstrate  a  method  for  interpreting  the  velocity  data  in  terms 
of  subsurface  volumetric  water  content. 

In  general,  ground  penetrating  radar  velocity  decreases  rapidly  with  depth  (Davis  and 
Annan,  1989).  This  is  primarily  a  result  of  increasing  water  saturation  with  depth.  Topp  et 
al.  (1980)  derived  an  empirical  relationship  between  dielectric  constant  of  soil  samples  and  the 
volumetric  water  content  of  the  samples,  where  volumetric  water  content  is  the  ratio  of  water 
volume  to  total  sample  volume.  Using  their  relationship,  or  similar  mixing  formulas,  we  can 
estimate  water  content  from  dielectric  constants  where  the  dielectric  constants  are  calculated 
from  the  interval  velocities  derived  from  normal  moveout  velocities.  This  interpretation  step 
increases  the  potential  uses  of  radar  profiling  in  ground  water  studies  and  contaminant  spill 
monitoring  [see  for  example  the  variations  in  radar  interval  velocity  caused  by  a  controlled 
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contaminant  spill  described  by  Brewster  and  Annan  (1994)].  Radar  interval  velocity  can  be 
used  to  estimate  water  content  when  the  subsurface  is  sufficiently  resistive  to  be  treated  as  a 
low-loss  media.  This  is  a  reasonable  assumption  where  radar  signals  penetrate  the  subsurface 
to  depths  of  meters  or  tens  of  meters.  In  partially  saturated  soils  with  constant  porosity, 
water  content  may  be  interpreted  as  an  indicator  of  saturation.  In  fully  saturated  soils, 
variations  in  water  content  can  be  interpreted  as  variations  in  porosity.  The  improvement  to 
the  radar  image  as  a  result  of  stacking  GPR  data  collected  in  the  CMP  geometry,  combined 
with  interpretation  of  the  subsurface  velocity  in  terms  of  water  content,  should  encourage 
the  development  of  GPR  systems  capable  of  acquiring  multi-offset  data  more  efficiently. 

2.  Multi-Offset  Data 

We  obtained  a  multi-offset  GPR  survey  from  Sensors  &  Software,  Inc.,  that  was  acquired  at 
the  Chalk  River  research  area,  operated  by  Atomic  Energy  of  Canada,  Ltd.  The  data  were 
acquired  in  cooperation  between  Sensors  &  Software,  Inc.,  the  Atomic  Energy  of  Canada, 
Ltd.,  and  the  UT-Dallas  Geophysical  Consortium.  The  data  were  collected  with  multiple 
source  antenna  to  receiver  antenna  offsets,  such  that  after  rearrangement,  1800  CMP’s  spaced 
every  0.25  m  were  defined.  The  data  were  acquired  using  a  pulse  EKKO  IV  digital  radar 
system  with  100  MHz  antennas.  A  detailed  description  of  the  acquisition  geometry  and  the 
data  recording  parameters  can  be  found  in  Fisher  et  al.  (1992).  For  these  data,  each  CMP 
gather  has,  on  average,  ten  traces  with  offsets  ranging  from  0.5  m  to  20.0  m.  The  standard 
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GPR  data  set  would  have  only  a  single  offset  trace  at  each  survey  station.  A  rule-of-thumb 
estimate  of  the  optimum  offset,  in  a  single  offset  survey,  is  that  the  antenna  separation  should 
be  about  20%  of  the  depth-of-interest  (Annan  and  Cosway,  1992).  In  Figure  1,  the  minimum 
(0.5  m)  offset  trace  profile  from  the  Chalk  River  data  set  is  displayed.  This  shows  the  data 
as  it  might  have  been  collected  in  a  single  offset  GPR  survey  at  this  site.  The  subsurface  at 
the  site  consists  of  bedrock  covered,  primarily,  by  aeolian  and  fluvial  sand  deposits  of  glacial 
origin  (Fisher  et  al.,  1992).  In  this  profile  we  see  a  strong  V-shaped  reflector  that  appears  to 
mark  the  base  of  the  sedimentary  section.  Also,  there  is  strong  loss  in  reflection  continuity 
between  about  100-400  ns  from  CMP  900  to  CMP  1800.  Improvements  to  the  radar  image 
in  these  areas,  in  particular,  are  observed  in  the  final  CMP  processed  data  profiles. 

Figure  2a  shows  a  number  of  CMP  gathers  with  traces  arranged  in  offset  order  from  0.5  m 
to  20.0  m  within  each  CMP.  Some  of  the  approximately  hyperbolic  arrivals  from  reflections 
are  clearly  distinguished,  but  the  signal-to-noise  is  low  in  this  raw  data.  Preprocessing 
steps  were  applied  to  the  data  to  prepare  it  for  velocity  analysis  and  stacking.  The  data 
were  bandpass  filtered  (20-200  MHz),  a  top  mute  (in  shot  gather  mode)  was  applied  to 
remove  the  direct  arrivals.  Amplitude  was  multiplied  by  the  square  of  sample  traveltime 
as  an  approximate  correction  for  geometric  spreading  loss.  An  f-k  dip  filter  was  applied  to 
common  shot  gathers  to  remove  some  reverse  dip  events  which  may  have  been  the  result  of 
scattering  from  off-axis  locations.  AGC  (automatic  gain  control)  scaling  is  used  for  display. 
Figure  2b  shows  the  same  CMP  gathers  after  these  processes  have  been  applied.  These 
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filtered  data  are  used  in  the  normal  moveout  velocity  analysis. 


3.  Normal  Moveout  Velocity  Estimation 

Since  the  data  in  this  survey  were  all  collected  with  the  CMP  geometry,  normal  moveout 
(NMO)  velocity  analysis  can  be  applied  at  any  or  all  of  the  CMP’s  to  define  the  subsurface 
velocity.  The  number  of  velocity  analyses  that  are  required  depends  on  how  strongly  the 
radar  velocity  varies  laterally.  There  is  no  rule  or  formula  for  determining  an  optimum 
number  of  velocity  analyses  for  processing  a  CMP  data  set,  so  this  parameter  must  be 
established  by  testing  the  data  itself.  In  this  example,  we  show  how  increasing  the  number 
of  velocity  analyses  affects  the  final  CMP  stacked  image. 

There  are  a  variety  of  schemes  used  in  NMO  velocity  analysis  (Yilmaz,  1987).  We  have 
used  a  semblance  amplitude  approach  where  the  data  in  the  CMP’s  are  normal  moveout 
corrected  and  stacked  using  a  range  of  trial  velocities.  The  amplitudes  versus  time  over 
the  whole  range  are  then  contoured  and  displayed  as  a  velocity  spectrum.  The  peaks  of 
the  amplitude  mapping  are  chosen  to  define  the  1-D  velocity  function  at  each  CMP  being 
analyzed.  After  a  number  of  velocity  functions  have  been  defined,  a  2-D  velocity  profile  is 
created  by  interpolation.  For  a  statistical  method,  like  semblance  mapping,  the  more  traces 
there  are  in  any  individual  CMP,  and  the  larger  the  offset  range  (within  the  small  spread 
approximation),  the  better  the  resolution  in  the  velocity  spectrum.  In  practice,  the  vertical 
resolution  of  NMO  velocity  analysis  depends  on  the  vertical  spatial  frequency  of  reflectors 
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strong  enough  to  produce  distinct  semblance  peaks.  Usually  we  pick  only  the  highest  ampli¬ 
tude  velocity  spectra  peaks  for  the  traveltime  versus  velocity  function.  Therefore,  the  NMO 
velocity  function  will  usually  have  fewer  velocity  layers  defined  than  there  are  reflections 
observed  in  the  data. 

When  the  number  of  traces  per  gather  is  small  and  somewhat  noisy,  we  expect  the 
velocity  spectra  to  have  a  relatively  poor  signal-to-noise  ratio.  In  order  to  improve  the 
picking  resolution  in  the  velocity  spectra  of  the  Chalk  River  data  we  combined  the  traces 
from  twenty  CMP  gathers  around  each  CMP  chosen  for  velocity  analysis.  A  typical  velocity 
spectrum  for  a  single  CMP  (no  combination)  is  shown  in  Figure  3a.  When  the  spectrum  of 
the  combined  20  CMP’s  is  found,  as  shown  in  Figure  3b,  the  velocity  peaks  are  much  better 
resolved  and  NMO  velocity  can  be  picked  with  more  confidence.  One  estimate  of  the  error 
in  the  velocity  pick  is  the  width  of  the  central  semblance  maximum.  For  example  at  550  ns 
the  peak  has  a  width  of  0.012  m/ns  which  gives  an  uncertainty  in  the  velocity  at  this  time 
of  about  ±8%.  To  an  extent,  resolution  in  the  spectra  seems  to  improve  with  increasing 
traveltime.  This  is  due  in  part  to  the  fact  that  for  near  surface  reflections,  the  further  offset 
traces  are  muted  and  do  not  contribute  to  the  analysis,  so  there  are  less  data  in  the  statistical 
analysis.  Secondly,  semblance  values  of  the  stack  change  more  slowly  as  velocity  increases. 
Since  radar  velocity  is,  in  general,  higher  in  the  near  surface,  we  expect  less  well-resolved 
peaks  in  the  spectra  of  near  surface  reflections. 

The  NMO  velocity  analyses  were  performed  at  regular  intervals  along  the  profile.  Figure 
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4  shows  the  interpolated  NMO  velocity  fields  as  the  number  of  velocity  analyses  is  increased. 
Figure  4a  is  the  simple  flat  layered  profile  that  results  when  only  one  CMP  (center)  is  an¬ 
alyzed.  Velocity  varies  from  about  .125  m/ns  near  the  surface  to  about  .06  m/ns  for  deep 
reflectors.  The  basement  and  the  overlying  sediment  interfaces  generate  strong  reflections 
which  are  laterally  continuous  across  the  image  and  allow  us  to  make  reliable  velocity  es¬ 
timates  in  the  sedimentary  section.  However,  the  reflections  observed  below  the  basement 
reflector,  probably  originating  from  either  fractures  or  foliations  (Fisher  et  al,  1992),  are  not 
sufliciently  strong  or  laterally  continuous  to  allow  accurate  velocity  estimates  to  be  made 
within  the  crystalline  basement.  The  deepest  reflectors  for  which  velocity  is  estimated  is  at 
about  650  ns  in  the  central  portion  of  the  data.  The  normal  moveout  velocity  determined 
for  the  last  strong  reflector  is  presumed  to  mark  the  top  of  the  crystalline  basement.  This 
velocity  is  used  to  do  the  NMO  correction  for  all  data  below  the  last  time  pick.  Figure 
4b  shows  the  2-D  velocity  profile  when  velocity  analyses  at  four  CMP’s  are  included.  It 
is  immediately  clear  that  there  is  significant  lateral  velocity  variation.  As  the  number  of 
velocity  analyses  increase  (Figures  4c-4e),  the  2-D  NMO  velocity  field  becomes  even  more 
heterogeneous,  although  the  general  background  layering  is  still  apparent.  Lateral  velocity 
variations  of  as  much  as  30%  are  encountered  in  a  traverse  across  the  more  detailed  velocity 
profiles. 
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4.  Stacked  Radar  Profiles 


Normal  moveout  corrections  calculated  from  the  NMO  velocity  field  remove  the  offset  de¬ 
pendence  of  the  reflection  traveltimes  such  that  the  data  can  be  treated  as  zero  offset  traces. 
After  the  correction  is  made,  the  traces  within  each  CMP  are  stacked  to  produce  a  single 
CMP  trace.  If  the  NMO  correction  is  done  with  the  correct  velocity  function,  the  stacked 
traces  have  an  improved  signal-to-noise  ratio. 

The  NMO  correction  is  based  on  the  assumption  that  the  subsurface  sampled  by  each 
CMP  can  be  adequately  modeled  as  a  sequence  of  horizontal  layers  with  uniform  interval 
velocity.  Steeply  dipping  reflectors  and  strong  velocity  gradients  test  the  applicability  of 
the  normal  moveout  and  stack  technique,  but  in  general  the  method  is  sufl&ciently  robust  to 
improve  data  quality  under  most  conditions.  If  lateral  velocity  variation  is  small  enough,  a 
single  NMO  velocity  function  can  be  used  to  calculate  the  moveout  correction  at  all  CMP’s. 
However,  where  lateral  velocity  variation  is  significant,  a  variable  NMO  velocity  field  defined 
by  functions  at  a  number  of  CMP’s  must  be  applied  to  obtain  the  best  results.  In  principle, 
an  NMO  velocity  defined  individually  by  velocity  analysis  for  every  CMP  should  yield  the 
most  accurate  result.  However,  a  spatial  limit  to  significant  (defined  by  the  interpreter) 
lateral  variation  in  NMO  velocity  is  usually  reached  at  some  multiple  CMP  spacing. 

NMO  corrections  were  computed  for  the  Chalk  River  data  using  each  of  the  velocity  fields 
in  Figure  4.  In  Figure  5a  we  show  the  stacked  profile  when  NMO  is  defined  by  the  single 
velocity  analysis  (Figure  4a).  Figure  5b  is  the  stacked  profile  result  when  all  35  velocity 
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analyses  (Figure  4e)  are  used.  First,  we  compare  these  to  the  standard  (no  stack)  profile 
shown  in  Figure  1.  The  stacked  data  profiles  show  a  strong  reduction  in  background  noise 
as  compared  to  single  offset.  Also,  some  deeper  reflections,  below  the  apparent  channel 
bottom  reflection,  have  been  greatly  enhanced,  in  particular,  the  reflector  at  about  625  ns 
from  CMP  250  to  CMP  1000.  Note  that  when  a  velocity  field  is  applied,  the  traveltime  to 
reflectors  changes  slightly.  This  is  a  consequence  of  the  normal  moveout  correction  which  is 
shifting  arrival  times  as  a  function  of  the  particular  velocity  field.  As  the  detail  in  the  applied 
velocity  field  is  increased,  the  structure  in  the  subsurface  image  becomes  more  accurate.  For 
example,  the  strong  reflector  at  600  ns  between  CMP  200  and  CMP  800,  becomes  more 
realistically  upwardly  concave  in  Figure  5b. 

Comparing  the  single  velocity  stack  (Figure  5a)  to  the  multiple  velocity  analyses  stack 
(Figure  5b),  we  observe  that  some  deeper  reflections  occur  down  to  750  ns  but  the  more 
clear  improvement  is  in  continuity  of  reflectors  throughout  the  section.  Overall,  stacking 
with  even  just  one  velocity  control  point  provides  the  majority  of  improvement  in  the  deep 
reflector  continuity  and  depth-of-imaging,  compared  to  no  stack.  With  increasing  detail  in 
the  velocity  field  the  primary  improvement  is  in  reflector  continuity  and  time  structure.  The 
fact  that  the  majority  of  the  increase  in  depth  of  imaged  reflectors  occurs  with  the  stack  from 
even  a  single  velocity  analysis  as  compared  to  no  stack,  suggests  that  it  is  the  multi-offset 
nature  of  the  stack,  in  particular  the  far  offsets,  that  provide  the  signal  from  the  deeper 
reflectors. 
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As  previously  noted,  the  data  between  about  200  to  400  ns  to  the  right  of  CMP  900 
shows  very  few  coherent  reflectors  in  the  standard  GPR  survey  shown  in  Figure  1.  In  fact, 
the  left  half  of  the  profile  is  quite  different  than  the  right  half  of  the  profile.  The  description 
of  the  survey  acquisition  by  Fisher  et  al.  (1992)  shows  that  there  is  a  bend  in  the  survey 
line  at  this  point.  However,  only  a  few  CMP’s  at  the  turning  point  are  affected  by  this 
change  in  geometry  and  does  not  explain  the  change  in  reflection  character.  With  the  NMO 
correction  and  stack,  the  region  on  the  right  changes  dramatically,  with  a  dramatic  increase 
in  the  number  of  coherent  reflectors.  When  only  one  velocity  is  applied  (Figure  5a),  this 
area  of  the  data  is  still  poorly  resolved.  However,  when  all  velocity  functions  are  included  in 
the  velocity  description,  the  stack  has  improved  continuity  of  a  number  of  reflections  from 
the  left  side  of  the  data  into  the  right  side.  Also,  note  that  the  deepest  reflection  on  the 
right  side,  now  appears  to  consist  of  a  series  of  step-like  events  carved  into  the  crystalline 
basement.  Although  not  shown  here,  only  small  differences  between  the  stack  based  on  18 
velocity  analyses  and  35  were  observed.  This  implies  that  lateral  velocity  resolution  for  these 
data  is  between  50  and  100  CMP’s  (12.5  to  25  m).  We  also  note  that  in  some  portions  of 
the  profiles  there  are  bursts  of  coherent  steeply  dipping  noise  that  may  be  further  evidence 
of  scattering  from  inhomogeneities  off  the  axis  of  the  profile.  These  could  be  removed  by 
applying  a  post-stack  dip  filter. 
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5.  Radar  Propagation  Approximations 


Radar  propagation  is  fundamentally  limited  by  the  conductivity  of  the  subsurface  medium. 
Only  in  low-loss  media  can  radar  signals  penetrate  deep  enough  to  provide  a  useful  subsurface 
image.  The  conductivity  of  common  soil  and  rock  materials  that  can  be  penetrated  to  useful 
survey  depths  is  less  than  10  mS/m  (Davis  and  Annan,  1989).  Useful  approximations  for 
describing  the  propagation  of  radar  signals  can  be  found  by  considering  the  time  harmonic 
solution  of  the  equation 


^  dE  d^E 

^  ^  a?  + 


(1) 


derived  from  Maxwell’s  equation  for  the  case  of  a  homogeneous  isotropic  medium.  E  is  the 
electric  field  amplitude  and  the  constants  of  proportionality,  e,  /a,  and  a  are  the  electric 
permittivity,  magnetic  permeability  and  conductivity  of  the  medium.  The  first  term  on  the 
right  side  of  the  wave  Eq  (1)  represents  conduction  of  charge  due  to  the  applied  electric  field. 
The  second  term  describes  the  displacement  of  charge  due  to  the  field. 

A  solution  to  Eq  (1)  of  the  form 


E  = 


yields  the  dispersion  relation 

P  =  iixou  — 

The  wavenumber,  k,  is  complex  and  can  be  written  as 

k  =  a  + i/3 


(2) 

(3) 

(4) 
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where  the  attenuation  constant 


a  =  (^\/l+tan2  5-lj  (5) 

and  phase  constant 

/?  =  Y  (^1  +  tan^  8  (6) 

are  both  real.  The  parameters,  c  and  /Cg  are  the  electromagnetic  velocity  in  free  space  and 
the  real  part  of  the  dielectric  constant  of  the  medium,  respectively.  The  loss  tangent  tan  6  is, 
for  a  plane  wave  propagating  with  the  angular  frequency  w,  the  ratio  of  conduction  current 
density  to  displacement  current  density. 

tan  S  =  — .  (7) 

uje 

Typically,  relative  susceptibility,  is  taken  to  be  unity  except  where  metallic  objects  or 
minerals  are  present  in  abundance  (Telford  et  at,  1976).  If  we  take  relative  susceptibility  to 
be  unity,  the  solution  to  the  equation  can  then  be  written  as 

E  =  (8) 

which  is  the  expression  for  a  damped  plane  wave  propagating  with  phase  velocity 

In  materials  that  are  good  conductors,  displacement  currents  are  negligible  compared 
to  conduction  currents,  and  Eq  (1)  reduces  to  a  diffusion  equation,  i.e.,  the  fields  do  not 
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propagate  as  electromagnetic  waves.  For  materials  with  low  conductivity,  and  when  the 
frequency  of  the  oscillating  electric  field  is  high  enough,  the  displacement  current  dominates 
over  the  conduction  current  and  electromagnetic  waves  will  propagate  (Stratton,  1941). 
Where  GPR  is  most  effective,  the  loss  tangent  is  very  small  and  the  diffusion  term  can 
usually  be  neglected.  In  this  case,  where  tan  6  <C  1,  the  phase  constant  reduces  to 

13  «  -^e  (10) 

c 


such  that 

c 


and  the  medium  attenuation,  a,  can  be  approximated  by 


a 

a  w  0.2  X  — 


(11) 


(12) 


where  o  is  in  units  of  mS/m.  At  a  depth  of  ^  meters  the  field  strength  is  reduced  to  ^  of  its 
original  amplitude.  Care  must  be  taken  in  interpreting  this  depth  in  terms  of  radar  propa¬ 
gation  distance.  Kong  (1990)  refers  to  this  as  the  penetration  depth  for  slightly  conducting 
media.  He  points  out  that  it  is  frequency  independent,  even  though  we  expect  penetration  to 
decrease  at  higher  frequencies,  and  is  effectively  a  volume  attenuation  due  to  wave  scattering 
within  the  medium.  It  is  tempting  to  refer  to  this  penetration  depth  as  ‘skin  depth’,  but  it  is 
preferable  to  reserve  this  nomenclature  for  the  frequency  dependent  penetration  depth  that 
is  found  from  the  attenuation  factor  in  highly  conducting  media.  Practical  definitions  of 
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depth  of  penetration  for  GPR  include  not  only  the  medium  attenuation  but  also  the  acqui¬ 
sition  system  transmission  and  receiving  characteristics  and  a  frequency  dependence  (e.g., 
Davis  and  Annan,  1989;  Daniels,  1989). 

This  discussion  has  centered  on  the  time  harmonic  solution  to  the  wave  equation  and 
as  such  has  not  considered  that  the  material  properties  are  also  frequency  dependent  and 
therefore  should  be  treated  as  complex  numbers.  In  the  frequency  range  of  GPR,  10-1000 
MHz,  the  frequency  dependence  of  the  dielectric  constant  and  conductivity  are  small,  in 
low-loss  media,  and  can  be  taken  to  be  constant  for  practical  purposes  (Davis  and  Annan, 
1989).  In  our  interpretation  we  treat  the  dielectric  constant  as  real  and  related  to  the  prop¬ 
agation  velocity  by  the  simple  expression  found  in  Eq  (11).  In  making  this  approximation 
we  recognize  that  frequency  dependence  as  well  as  many  other  factors  such  as  scattering 
loss,  source/receiver  antenna  power  and  transmission  characteristics,  and  ground  couphng 
are  not  being  accounted  for.  However,  our  assumption  is  that  these  factors  change  much 
more  slowly  than  dielectric  constant  if  reflections  are  observed  in  the  data. 

6.  Interval  Velocity  and  Water  Content 

To  interpret  the  NMO  velocity  field  derived  from  the  multi-offset  data,  it  is  necessary  to 
calculate  interval  velocities  and  find  the  relationship  of  radar  propagation  velocity  to  other 
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geoelectric  properties.  Interval  velocity  Vi^n  was  calculated  using  the  Dix  formula  (Dix,  1955). 


This  calculation  was  made  difficult  by  the  fact  that  radar  velocity  decreases  with  increas¬ 
ing  traveltime.  The  Dix  formulation  does  not  preclude  the  case  of  velocity  decreasing  with 
traveltime,  but  we  found  that  with  decreasing  velocity  the  numerator  inside  the  square  root 
of  (13)  can  be  negative  if  the  traveltime  interval  is  small  or  the  NMO  velocity  change  is  large. 
Where  this  situation  was  encountered,  we  used  the  average  interval  velocity  of  laterally  ad¬ 
jacent  intervals.  When  the  Dix  formula  is  applied  to  calculated  exact  NMO  traveltimes  and 
moveout  velocities  for  a  radar  velocity  model,  this  problem  is  not  observed.  This  implies 
that  the  Dix  inversion,  when  applied  to  radar  data,  is  very  sensitive  to  noise  in  the  velocity 
analysis.  In  consideration  of  this,  the  interval  velocity  was  subjected  to  smoothing  prior  to 
calculation  of  dielectric  constant  and  depth  values. 

The  final  step  in  our  interpretation  scheme  is  to  relate  the  calculated  radar  interval  veloc¬ 
ities  to  water  content.  Interval  velocity  derived  from  reflection  moveout  curves  is  the  group 
velocity.  As  discussed  in  the  previous  section,  if  we  assume  that  the  subsurface  properties  of 
the  media  penetrated  by  the  radar  signal  do  not  vary  much  in  the  limited  frequency  range 
of  the  radar  signal  being  used,  then  the  group  velocity  can  be  taken  to  be  approximately 
equal  to  the  phase  velocity.  Then  we  can  use  the  phase  velocity  Eq  (11)  to  convert  interval 
velocity  to  dielectric  constant. 

The  dielectric  constant  of  rocks  and  soils  is  very  sensitive  to  the  water  content  of  the 


16 


sample.  This  is  because  the  dielectric  constant  of  water  («„,  ~  80)  is  much  larger  than  the 
dielectric  constant  of  both  mineral  grains  (kj  ~  3  —  5)  and  air  («„  =  !)•  Therefore,  as 
the  water  content  of  a  rock  or  soil  sample  increases,  the  dielectric  constant  of  the  sample 
increases  and  its  velocity  decreases.  The  net  water  content  (6)  is  equal  to  the  product  of  the 
sample  porosity  (^)  and  the  saturation  (S.^) 

e  =  <j>S^.  (14) 

Therefore,  the  dielectric  constant  and  velocity  of  a  sample  is  determined  by  its  porosity  and 
saturation,  and  to  a  lesser  degree,  by  the  microgeometry  of  the  sample. 

Several  mixing  formulas  have  been  both  theoretically  and  empirically  developed  for  the 
dielectric  response  of  heterogeneous  mixtures  such  as  water  saturated  rocks  (see,  for  example, 
Jackson  and  O’Neill,  1986;  Dobson  et  ai,  1985).  One  such  mixing  formula  is  the  GRIM 
(complex  refractive  index  method)  equation,  which  is  often  used  in  the  interpretation  of 
electromagnetic  logging  data  (Schlumberger,  1991).  The  GRIM  equation  is  a  semi-empirical 
mixing  law  which  relates  the  dielectric  constant  of  the  sample  to  the  water  filled  porosity,  <f>. 
At  radar  frequencies,  where  tan  6  1,  the  GRIM  equation  can  be  expressed  as 

+ {1  -  (15) 

where  Ke,  Kw,  and  Kg  are  the  dielectric  constants  of  the  sample,  the  pore  water,  and  the  dry 
mineral  grains,  respectively  (Wharton  et  al,  1980).  Since  Kyj  and  Kg  can  be  taken  as  known 
constants,  the  GRIM  equation  can  be  used  to  estimate  the  porosity  of  water-filled  samples 
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from  the  measured  dielectric  constant  of  the  sample. 

In  a  similar  fashion,  the  Hanai-Bruggeman  (H-B)  mixing  formula  can  also  be  used  to 
estimate  the  porosity  of  a  water  saturated  rock  or  soil  sample  (see,  for  example.  Sen  et  ai, 
1981).  The  H-B  formula  is  obtained  by  a  theoretical  effective  medium  calculation  which 
expresses  the  complex  dielectric  response  of  the  sample  in  terms  of  the  complex  dielectric 
response  of  the  mineral  grains  and  the  pore  fluids  that  make  up  the  sample.  At  high  fre¬ 
quencies,  where  tan  5  <  1,  the  H-B  equation  yields  the  following  implicit  expression  for  the 
dielectric  constant  of  the  sample. 


/Cg  —  K/yjCf) 


(16) 


This  expression  can  be  explicitly  solved  for  the  sample  porosity  of  the  fully  saturated  sample. 


,K,n 


<I>HB  =  (-)’ 


(17) 


The  cementation  index,  m,  is  a  function  of  the  grain  shape  and  is  typically  found  to  vary 
between  1.5,  for  unconsolidated  sands  with  well-rounded  grains,  to  approximately  2.0  for 
well-cemented  sandstones  with  oblate  grain  shapes  (Jackson  et  al,  1978). 

Figure  6a  shows  the  dielectric  constant  of  a  water  saturated  sandstone,  predicted  by  the 
CRIM  and  H-B  miyirig  formulas,  as  a  function  of  the  sample  porosity.  In  these  simulations, 
we  have  taken  =  80,  Kg  =  4.5  and,  for  the  H-B  formula,  we  have  used  cementation 
exponents  of  m  =  1.5  and  tti  =  2.0.  For  all  three  curves,  the  dielectric  constant  increases 
with  increasing  sample  porosity.  At  a  fixed  porosity,  the  H-B  equation  predicts  a  lower 
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dielectric  constant  for  a  well-cemented  sandstone  (m  =  2.0)  than  for  a  poorly  cemented 
sandstone  with  round  grains  (m  =  1.5).  The  curve  corresponding  to  the  GRIM  formula, 
which  does  not  have  an  explicit  dependence  upon  the  sample  microgeometry,  is  bounded  by 
the  two  H-B  curves. 

Another  way  to  vary  the  water  content  in  a  rock  or  soil  sample  is  to  vary  the  level  of 
saturation.  In  this  analysis,  we  consider  a  partially  saturated  rock  as  a  three-phase  mixture 
consisting  of  mineral  grains,  pore  water,  and  air.  The  saturation,  Sy,,  is  the  fractional 
volume  of  the  pore  space  occupied  by  water,  and  therefore,  it  varies  from  Syj  =  0.0,  for 
an  unsaturated  sample,  to  Sy,  =  1.0,  for  a  fully  saturated  sample.  To  predict  variations  in 
the  dielectric  constant  of  partially  saturated  samples,  we  need  to  use  a  three-phase  mixing 
formula.  Both  the  GRIM  and  the  H-B  formulations  can  be  modified  to  predict  the  dielectric 
response  of  three-phase  mixtures. 

The  GRIM  formula  is  the  electrical  analog  to  the  Wiley  time-average  equation,  which  is 
frequently  used  to  predict  the  acoustic  velocities  of  porous  media  (see,  for  example,  Gueguen 
and  PaJciauskas,  1994).  In  these  time-averaged  equations  the  effective  slowness  (inverse  of 
velocity)  of  the  sample  is  obtained  by  adding  the  sample  phases  in  series  with  each  other. 
Therefore,  a  three-phase  mixture  of  grains,  water  and  air,  results  in  the  following  GRIM  type 
mixing  formula. 

\/ \/ (l  0) K>g  -j-  ^(l  ^ii;)\/ (^^) 

The  H-B  effective  medium  theory  can  also  be  used  to  predict  the  effective  dielectric  response 
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of  partially  saturated  rock  or  soil  samples  (Endres  and  Knight,  1992;  Samstag,  1992).  In  this 
analysis,  we  follow  the  procedure  of  Samstag  and  Morgan  (1991),  in  which  the  H-B  mixing 
formula  is  employed  twice;  first,  to  compute  the  effective  dielectric  response  of  the  water/ air 
mixture,  which  fills  the  pore  space,  and  second  to  compute  the  effective  dielectric  response 
of  the  total  rock,  which  is  obtained  by  mixing  the  mineral  grains  into  the  water/air  mixture. 

The  first  mixing  of  air  into  water  gives 

/  1  _  \ 

(19) 

\  «porc  / 


for  the  effective  dielectric  response  of  the  pore  ‘fluid’  mixture,  where  the  exponent,  mi,  is 
related  to  the  shape  of  the  air  bubbles.  The  second  mixing  of  the  mineral  grains  into  the 
water/air  mixture  gives 


Kq  —  l^pore^ 


(20) 


for  the  effective  dielectric  response  of  the  partially  saturated  sample,  where  m2  is  related 
to  the  shape  of  the  mineral  grains.  In  general,  the  geometric  factors,  mi  and  m2,  in  the 
first  and  second  embeddings,  may  be  different.  Furthermore,  mi  may  also  be  a  function  of 
saturation  (Endres  and  Knight,  1992;  Morgan  and  Samstag,  1991).  However,  in  this  analysis 
we  have  assumed  that  mi  is  independent  of  saturation,  and  equal  to  m2. 

Figure  6b  is  a  plot  of  the  H-B  and  CRIM  simulations  of  the  dielectric  constant  versus 
water  saturation  for  a  sample  with  40%  porosity.  As  in  Figure  6a,  the  H-B  simulations  with 
m  =  1.5  (spherical  inclusions)  and  m  =  2.0  (oblate  inclusions),  bound  the  dielectric  response 
obtained  from  the  three-phase  CRIM  mixing  formula. 
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In  practice,  it  is  impossible  to  derive  both  the  sample  porosity  and  water  saturation  from 
just  the  dielectric  constant  of  the  sample.  Topp  et  al.  (1980),  therefore,  used  a  wide  range 
of  soil  samples  with  varying  degrees  of  water  saturation  to  obtain  the  following  empirical 
relationship  between  the  measured  dielectric  response  of  the  samples  and  the  net  water 
content  {6  =  4)S.u,). 

Ke  =  3.03  +  9.3  ^  +  146.0  0^  -  76.7  e\  (21) 

We  note  that  in  the  Topp  et  al.  publication,  an  expression  is  also  given  for  finding  9  in  terms 
of  a  polynomial  in  Ke  (22). 

e  =  -5.3  X  10“2  +  2.92  X  10“^  Ke  -  5.5  x  10"'‘  kI  +  4.3  x  10~®  (22) 

In  comparing  crossplots  (not  shown  here)  of  Ke  and  9  calculated  from  the  two  expressions, 
we  found  that  for  the  range  of  water  content,  9  =  0.02  —  0.4,  the  expressions  are  in  good 
agreement,  however  at  ^  =  0.0  they  differ  by  about  one  dielectric  constant  unit. 

In  Figure  6c,  the  dielectric  response  of  the  Topp  equation  (21),  plotted  as  a  function  of 
water  content,  is  compared  to  the  H-B  responses  for  the  two-phase  and  three-phase  mixtures 
plotted  in  Figures  6a  and  6b,  respectively.  Recall  that  in  the  two-phase  H-B  simulations, 
the  samples  are  fully  saturated  and  the  porosity  is  varied  from  0  —  40%  {9  =  cp),  whereas,  in 
the  three-phase  H-B  simulations,  porosity  is  held  constant  at  0  =  0:'4  and  Sju  is  varied  from 
0.0  to  1.0  {9  =  0.4  S^). 

Note  that  at  high  water  content,  both  the  two-phase  and  the  three-phase  H-B  simulations 
give  comparable  results,  with  the  m  =  1.5  case  giving  larger  dielectric  constants  than  the 
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m  =  2.0  case.  At  low  water  content,  the  two-phase  H-B  simulations  converge  to  Kg  —  2.5 
which  corresponds  to  the  dielectric  constant  of  an  unsaturated  sandstone. 

The  Topp  equation  was  obtained  by  fitting  dielectric  measurements  made  for  a  wide 
variety  of  soil  types  with  varying  porosities,  water  saturations,  and  cementation  exponents. 
Therefore,  it  is  not  possible  to  make  a  direct  comparison  between  the  response  of  the  Topp 
equation  and  the  responses  of  the  H-B  model.  However,  it  is  interesting  to  note  that  at  low 
water  content,  the  Topp  equation  is  in  good  agreement  with  the  three-phase  H-B  model, 
with  m  =  2.0,  and  at  high  water  content  the  Topp  equation  is  in  good  agreement  with  the 
two-phase  and  three-phase  H-B  models,  with  m  =  1.5.  One  possible  explanation  for  this 
may  be  that  in  partially  saturated  samples,  the  microgeometry  of  the  water/air  mixture, 
that  fills  the  pore  space,  may  actually  vary  with  the  level  of  water  saturation  (Endres  and 
Knight,  1992).  The  Topp  curve  may  actually  indicate  that  the  cementation  exponent  for  the 
water/air  mixture  in  partially  saturated  samples  decreases  with  increasing  saturation. 

In  the  following  analysis,  we  use  the  Topp  equation  to  transform  the  dielectric  constants, 
obtained  from  the  interval  velocities,  into  estimates  of  the  subsurface  water  content.  As 
the  curves  in  Figure  6c  imply,  the  transformation  of  dielectric  constant  into  water  content 
is  a  function  of  the  rock  porosity,  saturation  and  microgeometry,  which  are  all  unknown. 
Therefore,  the  transformation  will  not  be  exact,  but  the  Topp  equation  should  adequately 
reflect  the  relative  changes  in  the  subsurface  water  content.  Note  that,  if  the  GPR  data  were 
collected  at  multiple  frequencies,  it  might  be  better  to  use  the  Hanai-Bruggeman  equation 
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which,  in  its  complete  form,  is  frequency  dependent.  It  should  also  be  noted  that  the  actual 
measured  property,  interval  velocity,  changes  more  slowly  with  water  content  as  velocity 
decreases.  Since  radar  velocity  in  general  decreases  with  depth,  decreasing  sensitivity  of  the 
inversion  to  variations  in  water  content  is  expected  to  occur  at  greater  depths. 

After  calculating  the  dielectric  constants  from  the  smoothed  interval  velocities  using 
Eq  (11),  we  use  the  Topp  Eq  (22)  to  estimate  the  subsurface  water  content  shown  in  Figure 
7.  We  also  use  the  interval  velocities  to  convert  the  profile  from  time  to  depth.  Before  making 
this  conversion,  the  velocity  estimates  were  binned  into  20  ns  intervals  and  then  smoothed 
with  a  nine-point  averaging  window.  The  resulting  water  content  depth  section  in  Figure  7 
clearly  shows  a  zone  of  high  water  content  in  the  central  portion  of  the  image.  The  bottom 
of  the  depth  section  reflects  the  irregular  subsurface  topography  of  the  crystalline  basement. 
Below  the  basement  reflector,  we  were  not  able  to  estimate  velocity  and  consequently  could 
not  estimate  water  content  within  the  crystalline  basement. 

The  hydrogeology  of  the  Twin  Lake  Site  at  the  Chalk  River  Nuclear  Laboratories,  where 
these  data  were  acquired,  has  been  extensively  characterized  with  a  large  network  of  boreholes 
and  several  hydrogeological  and  geophysical  studies  (Thomas,  1989;  Catto  et  al,  1989;  Kiley 
and  Annan,  1989).  Fisher  et  al.  (1992)  used  stratigraphic  information  obtained  from  five 
boreholes,  located  in  the  vicinity  of  the  radar  survey,  to  interpret  the  processed  radar  data. 
An  adapted  version  of  their  interpretation  is  shown  in  the  lower  part  of  Figure  7.  The  top 
stratigraphic  layer  is  composed  of  unconsolidated  aeolian  sands  predominantly  medium-fine 
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in  grain  size  with  lenses  of  fine  sand  and  silt  stringers  (indicated  by  the  dotted  lines  in  the 
figure).  The  bottom  of  the  aeolian  deposits  (indicated  by  the  dark  solid  line  in  the  figure) 
is  marked  by  a  thin  (0.5  m  to  2.0  m)  sand  zone  which  is  enriched  in  garnet  and  other  heavy 
minerals.  The  sand  layers  below  the  garnetiferous  sand  are  fluvial  deposits,  consisting  of 
unconsolidated  sand  of  fine  to  medium  grain  size,  with  a  silty  clay  layer  located  toward  the 
bottom  left  side  of  the  sedimentary  section.  The  unit  below  the  silty  clay  layer  was  labeled 
as  ‘unknown’  by  Fisher  et  al.  because,  although  it  was  observed  in  the  radar,  it  was  not 
penetrated  by  a  borehole. 

The  radar  survey  was  run  along  a  road,  which  at  CMP  805  changed  direction  and  began 
to  climb  up  a  slight  gradient.  There  is  approximately  10  m  of  elevation  gain  between  CMP 
805  and  CMP  1800.  Topographic  corrections  have  not  been  made  to  either  the  water  content 
depth  section  or  the  geological  cross-section.  If  topographic  corrections  were  made  to  the 
water  content  section,  the  zone  of  high  water  content  between  CMP  805  and  CMP  1800 
would  be  pulled  up  by  the  elevation  change.  In  this  case,  the  high  water  content  would  show 
a  downdip  trend  from  right  to  left  in  the  profile.  This  direction  corresponds  to  the  direction 
of  regional  groundwater  flow  (Killey  and  Annan,  1989).  However,  it  is  also  apparent  that 
there  will  still  be  some  isolation  of  the  high  water  content  zone  in  the  center  of  the  cross- 
section  that  may  be  indicative  of  water  table  mounding  where  the  profile,  which  runs  from 
northwest  on  the  left  to  southeast  on  the  right,  crosses  between  Twin  Lake  to  the  south  and 
a  wetland  area  to  the  north. 
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The  unknown  stratigraphic  unit,  identified  by  Fisher  et  al.  in  the  radar  profile,  seems  to 
have  a  higher  water  content  than  the  overlying  silty-clay  unit.  Since  it  has  water  content 
similar  to  the  other  fluvial  sands  in  the  section,  it  may  also  be  a  fluvial  sand.  Our  water 
content  depth  section  indicates  a  maximum  depth  to  basement  of  approximately  30  m, 
whereas  Fisher  et  al.  indicate  a  maximum  depth  to  basement  of  about  25  m.  Killey  and 
Annan  (1989)  report  basement  depths  in  this  area  ranging  from  0  to  30  meters,  making  it 
impossible  to  be  certain  which  depth  calculation  is  correct.  The  discrepancy  in  basement 
depth  between  our  analysis  and  that  of  Fisher  et  al.  corresponds  to  about  a  15%  difference  in 
average  velocity  to  basement.  We  suspect  that  much  of  the  error  is  in  the  interval  velocities 
of  the  shallower  portions  of  the  section,  where  NMO  velocity  picking  was  most  diflflcult.  In 
any  case,  it  is  advisable  to  consider  the  water  content  values  as  showing  relative  changes 
rather  than  relying  strongly  on  the  absolute  values. 

7.  Conclusions 

Ground  penetrating  radar  surveys  collected  with  the  CMP  multi-offset  geometry  yield  im¬ 
proved  subsurface  images  over  single  offset  surveys.  We  have  shown  that  the  CMP  profile 
is  improved  as  the  number  of  velocity  analyses  is  increased.  This  leads  to  the  conclusion 
that  lateral  variation  in  the  radar  propagation  velocity  profile  can  be  significant  even  over 
the  limited  range  of  a  typical  GPR  survey.  The  CMP  stacking  process  yields  improved 
depth-of-imaging  compared  to  a  single  offset  survey  while  detailed  velocity  analysis  yields 
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improvement  to  continuity  of  reflectors  throughout  the  stacked  proflle. 

In  this  study,  we  have  attempted  to  connect  the  practical  measurement  of  radar  velocity 
from  the  multi-offset  data  with  the  theoretical  and  laboratory  relationships  between  water 
content  and  dielectric  constant.  We  have  shown  a  practical  approach  to  estimating  water 
content  from  these  velocities  based  on  the  assumption  that  in  low-loss  media,  radar  prop¬ 
agation  velocity  can  be  expressed  as  a  function  of  the  real  part  of  the  dielectric  constant. 
We  applied  this  method  to  the  Chalk  River  multi-offset  GPR  data  and  calculated  a  water 
content  cross-section  in  depth.  Comparison  of  this  result  to  published  geologic  and  hydroge¬ 
ologic  descriptions  of  the  survey  area  suggests  that  our  water  content  depth  section  is  indeed 
reasonable.  Based  on  similarity  of  water  contents,  we  suggest  that  the  sedimentary  unit 
described  by  Fisher  et  al.  (1992)  as  ‘unknown’  is,  in  fact,  a  fluvial  sand  deposit. 

Improvement  to  this  approach  could  be  made  by  increasing  the  number  of  source  to 
receiver  offsets  recorded  for  each  CMP  which  would  improve  the  NMO  velocity  analysis. 
Also,  we  suggest  that  a  test  of  this  interpretation  method  should  be  done  by  measuring 
subsurface  water  content  in  an  area  while  acquiring  a  multi-offset  radar  survey. 
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Figure  2a:  Radar  CMP  gathers  in  offset  order  before  preprocessing.  Offsets  vary  from  0.5 
m  to  20.0  m  within  each  CMP. 
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Figure  2b:  Radar  CMP  gathers  after  preprocessing.  Offsets  vary  from  0.5  m  to  20.0  m  within 
each  CMP. 
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Figure  3:  (a)  Velocity  spectrum  of  CMP  755.  (b)  Velocity  spectrum  of  the  combination  of 
20  CMP’s  centered  on  CMP  755. 


Figure  5a;  CMP  stacked  radar  reflection  profile  using  the  velocity  profile  from  1  CMP, 
corresponding  to  Figure  4a. 
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Figure  5b:  CMP  stacked  radar  reflection  profile  using  velocity  profile  from  35  CMP’s,  cor¬ 
responding  to  Figure  4e. 


POROSITY 

Figure  6a:  Dielectric  constant  as  a  function  of  porosity  (0)  as  predicted  by  the  GRIM  (-) 
and  Hanai-Bruggeman  (...)  two-phase  mixing  laws.  The  cementation  exponents  (m)  for 
the  H-B  simulations  are  indicated  on  the  figure. 


Sw 


Figure  6b:  Dielectric  constant  as  a  function  of  water  saturation  as  predicted  by  the 
GRIM  (-)  and  Hanai-Bruggeman  (-  -)  three-phase  mixing  laws.  The  sample  porosity  is 
held  constant  6  =  0.4  and  the  cementation  exponents  (m)  for  the  H-B  simulations  are 
indicated  on  the  figure. 
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Figure  6c:  Dielectric  constant  as  a  function  of  water  water  {9  =  (j)Sw)  as  predicted  by  the 
Hanai-Bruggeman  effective  medium  theory  for  both  two-phase  (...)  and  three-phase  (-  -) 
mixtures,  and  the  empirical  Topp  et  al.  (1980)  equation  (-).  The  cementation  exponents 
(m)  for  the  H-B  simulations  are  indicated  on  the  figure. 
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Figure  7:  Comparison  of  the  water  content  depth  section  (middle),  estimated  from  radar 
interval  velocities,  to  the  geologic  interpretation  of  Fisher  et  al  (1992).  An  approximate 
topography  profile  is  shown  at  the  top  of  the  figure.  In  the  geologic  cross-section,  the 
dotted  lines  indicate  thin  silt  layers,  and  the  dark  solid  line  indicates  a  garnet  sand  which 
appears  to  define  the  contact  between  the  aeolian  sands  near  the  surface  and  the  deeper 
fluvial  deposits. 
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HIGH-RESOLUTION  SHALLOW  SEISMIC  STRUCTURE 
IMAGING  USING  GRID-BASED  NONLINEAR  REFRACTION 

TRAVELTIME  TOMOGRAPHY 


Summary 

We  develop  a  nonlinear  refraction  traveltime  tomography  technique  that  can  rapidly  recon¬ 
struct  seismic  velocity  distribution  in  a  large  2-D  gridded  model.  For  the  traveltime  and 
raypath  calculations,  we  improve  a  graph-theoretical  method  by  placing  nodes  on  each  grid 
boundary  with  an  optimized  distribution  pattern.  This  method  greatly  reduces  errors  in 
traveltime  and  raypath  calculations  for  complex  velocity  models.  We  pose  a  nonlinear  objec¬ 
tive  function  for  inversion,  which  includes  not  only  the  traveltime  misfit  and  model  curvature 
roughness  norm,  but  also  traveltime  gradient  misfit  norm.  Our  numerical  experiments  show 
that  fitting  the  gradients  of  the  traveltime  curves  in  addition  to  the  traveltime  itself  can  bet¬ 
ter  resolve  the  velocity  contrasts  across  interfaces.  We  apply  Newton’s  method  to  minimize 
the  nonlinear  objective  function  iteratively.  The  use  of  Tikhonov  regularization  allows  to 
perform  a  global  inversion  and  reconstruct  the  whole  model  with  minimum  model  curvature 
roughness. 

We  apply  this  refraction  traveltime  tomography  technique  to  image  shallow  bedrock  to¬ 
pography  at  a  coastal  site  near  Boston,  Massachusetts.  We  will  demonstrate  that  performing 
nonlinear  tomography  can  effectively  reconstruct  the  subsurface  image  for  complex  velocity 
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structure. 


1.  Introduction 

The  seismic  refraction  method  is  a  very  useful  tool  for  investigating  problems  at  different 
scales,  from  environmental  geophysics  and  oil  exploration  to  crustal  structure  studies.  Con¬ 
ventionally,  seismic  refraction  data  are  acquired  with  “forward”  and  “reverse”  geometry,  and 
the  interpretation  involves  reciprocity  and  some  ideas  similar  to  those  in  reflection  migra¬ 
tion.  These  include  generalized  reciprocal  method  (Palmer,  1980),  wavefront  reconstruction 
method  (Aldrige  and  Oldenburg,  1992),  and  wavefleld  extrapolation  method  (Hill,  1987). 
However,  all  these  methods  can  only  be  applied  to  areas  where  seismic  velocity  structures 
are  relatively  simple.  Tomography  methods,  using  an  approach  to  calculate  traveltime  and 
raypaths  in  a  gridded  model  and  applying  an  inversion  technique  to  reconstruct  seismic  ve¬ 
locities,  have  been  extensively  studied  and  developed  for  crosshole  geometry,  but  little  has 
been  done  for  refraction  problems.  Docherty  (1992)  demonstrates  a  2-D  refraction  travel¬ 
time  tomography  that  applies  to  a  two-layer  model  for  imaging  shallow  weathering  layer. 
Hole  (1992)  iteratively  solves  a  linearized  backprojection  problem  with  somewhat  ad  hoc 
smoothing  constraints  for  areas  in  the  model  that  do  not  contain  any  raypath.  Ammon  and 
Vidale  (1993)  conduct  refraction  traveltime  tomography  using  a  finite-difference  method 
(Vidale,  1988)  for  traveltime  calculations  without  rays  involved,  but  have  to  repeat  forward 
calculations  for  numerous  times  in  order  to  obtain  sensitivity  information.  Qin  et  al.  (1993) 
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improved  finite-difference  method  for  better  accuracy  but  with  a  higher  cost  of  computation 
time.  To  avoid  local  minima  in  inversion,  they  apply  an  ad  hoc  horizontal  smoothing  opera¬ 
tor  in  the  model  space  and  also  downward  extrapolate  the  gradient  of  the  traveltime  misfit 
function. 

In  this  report,  we  introduce  a  different  approach  for  rapidly  and  accurately  conducting 
nonlinear  grid-based  refraction  traveltime  tomography.  We  account  for  the  nonlinearity  of 
the  traveltime  inversion  problem  by  updating  raypaths  for  perturbed  model.  Our  forward 
traveltime  and  raypath  calculations  are  performed  using  a  graph-theoretical  method  based 
on  Saito  (1989)  and  Moser  (1990)  but  with  significant  improvement  in  accuracy  and  effi¬ 
ciency.  We  pose  a  nonlinear  objective  function  that  is  particularly  meaningful  for  refraction 
traveltime  problems,  i.e.,  including  the  misfit  norm  of  the  gradients  of  traveltime  curves, 
which  contain  independent  slowness  information  in  addition  to  the  traveltime.  As  opposed 
to  applying  ad  hoc  model  smoothing  to  keep  the  inversion  stable,  we  solve  the  nonlinear 
traveltime  tomography  problem  using  the  regularization  method  of  Tikhonov.  Tikhonov 
regularization  takes  a  more  explicit  approach  by  damping  spatial  curvature  of  the  model 
function  and  is  applied  to  perform  a  global  inversion  in  the  sense  of  reconstructing  the  whole 
model.  The  inverse  problem  is  solved  by  Newton’s  method  (Tarantola,  1987),  and  conjugate 
gradient  approach  (Hestenes  and  Stiefel,  1952)  is  applied  to  rapidly  update  model  changes 
in  each  inversion  iteration.  We  demonstrate  in  real  cases  that  this  approach  can  produce 
a  stable,  unique  solution  of  the  nonlinear  refraction  traveltime  problem.  For  a  model  with 


41 


250  X  100  grids,  and  12  sources  and  24  receivers  applied,  tomography  can  be  solved  in  about 
7  min  on  a  DEC  alpha  workstation  starting  from  an  arbitrary  homogeneous  model. 

2.  Traveltime  and  Raypath  Calculations 

For  traveltime  inversion,  a  forward  modehng  procedure  is  needed  not  only  for  the  traveltime 
calculation  but  also  for  the  raypath  determination.  Raypath  tells  model  sensitivity  to  the 
traveltime  measurements.  One  can  use  the  finite-difference  method  to  solve  eikonal  equation 
for  traveltime  calculation  (Vidale,  1988;  Qin  ct  a,l.,  1993),  but  additional  efforts  are  needed 
to  find  raypaths.  We  chose  a  graph-theoretical  method  that  can  calculate  traveltime  and 
determine  raypaths  simultaneously  (Saito,  1989;  Moser,  1990). 

In  the  graph-theoretical  method,  nodes,  where  traveltime  and  raypaths  are  tracked,  are 
often  placed  at  the  corners  of  the  model  grids  (Saito,  1989;  Moser,  1990;  Mandal,  1992; 
Matarese,  1993)  or  regularly  on  the  grid  boundaries  (Nakanishi  and  Yamaguchi,  1986;  Moser, 
1990).  We  chose  the  one  which  has  nodes  distributed  on  the  grid  boundaries  and  straight 
rays  do  not  cross  grid  interfaces  unless  the  local  medium  is  homogeneous.  However,  instead 
of  equally  placing  nodes  on  the  grid  boundaries,  we  optimize  node  distribution  in  the  graph 
template  so  that  the  angle  propagation  error  is  greatly  reduced.  This  is  particularly  im¬ 
portant  for  refraction  traveltime  calculations  since  head  waves  are  generated  precisely  at 
critical  angles.  For  a  discretized  model,  there  are  two  independent  errors  involved  in  graph 
method;  error  due  to  finite  grid  size  and  error  due  to  propagation  angle  sampling  (Moser, 
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1990;  Klimes  and  Kvasnicka,  1993).  Decreasing  grid  size  cannot  reduce  propagation  angle 
error,  and  the  angle  error  is  associated  with  the  largest  angle  diflFerence  that  can  be  sampled. 
Therefore,  minimizing  angle  differences  for  a  given  number  of  nodes  simply  reduces  angle 
error.  Figure  1  illustrates  a  graph  template  that  contains  two  nodes  on  each  grid  boundary. 
The  node  distribution  is  optimized,  and  all  the  propagation  angles  are  equal  (22.5°).  If  the 
number  of  nodes  on  each  grid  boundary  is  more  than  two,  then  the  optimized  propagation 
angles  cannot  be  exactly  equal,  but  have  an  overall  minimum  difference. 

In  the  graph  method,  the  most  time  consuming  job  is  the  determination  of  shortest 
traveltime  among  some  selected  nodes.  We  use  an  INTERVAL  method  which  is  about  four 
to  five  times  faster  than  the  heap  sorting  algorithm  (Klimes  and  Kvasnicka,  1993). 


3.  Inversion  Method 


We  pose  a  nonlinear  objective  function  which  is  physically  meaningful  for  the  refraction 
traveltime  tomography: 


$(m)  =  (l-u;)||Q(d-G(m))|p+a;||r>,(d-G(m))||'  +  r||Lm|l2 
=  (1  -  uj)\\md  -  ffiG{m)\f  +  Lo\\md  -  7hG{m)\f  +  T\\Lm\f, 


(1) 
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def 

d 

rrid 

l{m) 
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dd 
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dx 
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where  d  is  the  trciveltinie  datsij  is  the  calculated  traveltime  data  for  current  model  tn. 

Cd  is  an  operator  that  divides  a  traveltime  by  the  corresponding  ray  length  I  and  returns 
an  average  slowness  frid  along  the  raypath.  Dx  is  a  differential  operator  for  traveltime 
with  respect  to  the  surface  geophone  spacing,  and  rhd  =  Dxd  returns  the  gradients  of  the 
traveltime-distance  curves,  L  is  a  Laplacian  operator  for  performing  Tikhonov  regularization, 
r  is  a  smoothing  trade-off  parameter,  and  a;  is  a  weighting  factor  between  data  misfit  norm 
and  data  gradient  misfit  norm.  The  inclusion  of  the  second  term  in  the  objective  function 
is  particularly  meaningful  for  refraction  traveltime  problem.  It  gives  an  h  norm  of  the 
apparent  slowness  misfit,  which  is  independent  to  the  information  from  traveltime  itself.  In 
fact,  picking  reliable  first  arrivals  from  seismograms  has  to  be  based  on  not  only  the  first 
break  on  the  individual  trace  but  also  on  the  systematic  moveout  of  a  traveltime  curve  across 
all  the  traces.  This  is  because  individual  trace  may  be  contaminated  by  incoherent  noise  or 
significant  changes  in  seismic  wavelengths  may  occur  due  to  site  effects  and  attenuation.  As 
one  can  see,  all  three  terms  as  defined  in  the  objective  function  (2)  are  explicitly  related  to 
model  slowness  in  one  way  or  another.  The  first  term  gives  misfit  of  the  average  slowness 
between  sources  and  receivers  along  raypaths;  the  second  term  represents  misfit  of  apparent 
slowness  along  the  surface;  and  the  third  term  calculates  the  curvature  roughness  of  the 
model  slowness. 

We  minimize  the  nonlinear  objective  function  (1)  using  Newton’s  method  with  some 
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modification: 


{{l-u;)AlAk-\-u}BlBk+TL'^  L+ekI)l^mk  =  {l-u})Al{md-mG{rnk))AujBl{md-mG{rnk)))-TL^  Lruk, 


(5) 

A  ^ 

drriG  _ 

1  dG 

(6) 

dm 

l{m)  dm' 

B 

drhc  _ 

d'^G 

(7) 

dm 

dmdx' 

rrik+i  =  mk  +  Aruk, 

k  =  1,2, 3, ...,  N, 

(8) 

where  Ak  contains  raypath  lengths  across  each  model  grid,  Bk  contains  the  difference  of 
raypath  lengths  in  each  model  grid  associated  with  the  traveltime  picks  at  two  adjacent 
geophones.  We  found  that  refraction  traveltime  inversion  behaves  similar  to  many  other 
nonhnear  geophysical  inversion  problems  in  the  sense  that  large  nonlinearity  occurs  in  the 
early  inversion  stage  due  to  a  poor  starting  model,  and  it  approaches  linearized  when  the 
model  is  close  to  the  true  solution.  Therefore,  we  include  a  variable  damping  term  tkl  in 
the  left-hand  side  in  Eq  (5),  and  define  Cfc  =  a  x  rhs,  where  a  is  an  empirical  parameter 
(about  0.01)  and  rhs  is  the  rms  misfit  norm  of  the  right-hand  side  in  Eq  (5).  If  the  objective 
function  is  not  minimized  well  and  quite  nonlinear,  then  rhs  is  large,  and  a  large  damping 
is  applied  and  small  model  updates  are  allowed.  With  the  inversion- proceeding  further  and 
rhs  decreasing,  a  smaller  €k  drives  the  convergence  speed  faster  in  the  latter  stage. 

We  apply  Tikhonov  regularization  in  our  inversion  as  given  by  the  third  term  in  Eq  (1) 
or  (2),  and  L  is  defined  as  the  second-order  difference  operator.  For  an  allowed  data  misfit 
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norm  level,  there  should  be  one  unique  trade-off  parameter  r.  Therefore,  experiments  with 
different  r  should  be  performed  in  order  to  find  the  appropriate  one.  This  selection  will  be 
based  on  the  convergence  curve  of  the  data  misfit  norm  versus  iterations. 

4.  Field  Data 

At  a  coastal  site  near  Boston,  Massachusetts,  a  small-scale  seismic  refraction  survey  was 
conducted  to  map  bedrock  topography.  The  goal  of  this  survey  was  to  locate  those  areas 
where  bedrock  is  deep  so  that  construction  of  a  new  storm-drainage  system  may  proceed 
without  costly  blasting.  The  environment  at  the  working  site  was  quite  unusual  as  the  survey 
area  was  covered  by  sea  water  during  high-tide  period,  and  exposed  only  for  about  1  or  2 
hours  during  the  low  tide  each  day.  Further  details  of  the  project  are  reported  by  Kutrubes 
et  al.  (1996). 

Figure  2  shows  the  seismic  waveforms  recorded  from  a  forward  shot  and  a  reverse  shot 
on  line  1.  Data  from  both  shots  show  relatively  delayed  first  arrivals  between  receivers  7  and 
14,  although  the  topography  along  line  1  is  flat.  Moreover,  the  amplitudes  of  these  delayed 
first  arrivals  are  relatively  small.  The  evidence  simply  suggests  that  a  low-velocity  zone 
with  strong  seismic  attenuation  exists  beneath  receivers  7  to  14.  It  becomes  more  obvious 
when  we  placed  sources  (using  hammer  and  air  gun)  at  locations  between  the  receiver  -7 
and  14,  such  effects  were  then  observed  at  all  the  traces.  Figure  3a  displays  traveltime  data 
from  12  shots  along  line  1.  Corrections  for  trigger  time  were  made  for  a  few  shots  based  on 
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reciprocity.  These  traveltime  curves  do  not  suggest  simple  velocity  structure,  rather,  indicate 
complexity  in  the  shallow  seismic  media.  In  another  case.  Figure  4  displays  waveform  data 
recorded  from  survey  line  2  and  Figure  3b  shows  traveltime  data  from  12  shots.  In  contrast 
to  line  1,  this  is  a  case  that  demonstrates  the  influences  due  to  high-velocity  anomalies  in 
the  shallow  structure,  i.e.,  an  intermediate  velocity  zone  seats  on  the  bedrock  and  outcrops 
the  surface  in  the  central  area.  Using  these  traveltime  data,  we  performed  tomography 
studies  with  models  consisting  of  250  x  100  grids.  Our  results  are  presented  in  Figure  5. 
The  calculated  traveltime  data  corresponding  to  our  final  solutions  are  plotted  in  Figure  2 
(grey  dots).  As  one  can  see  from  these  cases,  although  recorded  data  show  complexities, 
the  resolved  bedrock  topography  is  quite  simple.  The  complexities  are  mostly  due  to  the 
shallow  velocity  structures.  These  studies  simply  justify  the  importance  of  using  nonlinear 
tomography  approach  for  imaging  the  shallow  earth. 

5.  Conclusions 

In  this  report,  we  described  an  accurate  and  efficient  approach  for  calculating  traveltime 
and  raypaths  for  any  velocity  models.  We  introduced  a  nonlinear  traveltime  tomography 
method  that  uses  first  arrivals  and  resolves  velocities  in  a  gridded  model.  Using  Tikhonov 
regularization,  we  obtain  a  global  solution  by  nonlinear  inversion.  Finally,  the  validity  of  our 
approach  was  proven  by  applications  to  real  data. 
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Optimized  Node  Distribution 


a1  =  a3  =  0.29289dx 

a2  =  0.41421dx 
0  =  22.5° 


Figure  1:  An  optimized  graph  template  in  which  the  differences  of  ray  angles  are  minimized. 
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Forward  Shot 


Reverse  Shot 
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Figure  2:  Waveform  data  from  forward  and  reverse  shots  on  line  1.  The  influences  due  to  a 
shallow  low  velocity  zone  are  shown  between  trace  7  and  14.  Time  interval  is  5  ms. 


a.  Calculated  Traveltime  and  Field  Data  for  Line  1 
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b. Calculated  Traveltime  and  Field  Data  for  Line  2 
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Figure  3:  Field  traveltime  data  (curves)  from  line  1  (a)  and  line  2  (b)  and  the  calculated 
traveltimes  (grey  dots)  for  the  resolved  models  shown  in  Figure  5. 
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Figure  5:  Tomographic  imaging  results  for  refraction  survey  line  1  and  line  2. 
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3-D  D.C.  ELECTRICAL  RESISTIVITY  INVERSION  WITH 
APPLICATIONS  TO  IMAGE  THE  ABERJONA  WATERSHED 

Summary 

The  task  of  d.c.  electrical  resistivity  inversion  is  to  solve  for  the  subsurface  resistivity  dis¬ 
tribution  from  the  measurements  of  electrical  potential  on  the  surface  of  the  earth  or  in 
boreholes.  This  problem  is  ill-posed  in  the  sense  that  an  infinite  number  of  solutions  to  the 
inverse  problem  exist  due  to  incomplete  and  uncertain  data.  To  overcome  the  ill-posedness 
solutions  must  incorporate  prior  information  or  a  smoothness  constraint  in  order  to  be  phys¬ 
ically  meaningful  and  stable  against  noise  in  the  data.  In  this  paper  we  investigate  inversion 
algorithms  based  on  Tikhonov’s  regulajization  method,  which  solves  a  minimization  problem 
to  find  models  that  fit  the  data  and  also  have  minimum  structure.  We  investigate  different 
smoothness  constraints  and  find  that  minimization  of  second-order  spatial  derivatives  of  the 
resistivity  to  be  most  effective,  particularly  in  avoiding  unstable  features  near  the  source  and 
receiver  electrode  locations. 

To  implement  the  regularized  inversion,  we  developed  a  numerical  algorithm  for  miiiimiz- 
ing  the  objective  functional  based  on  nonlinear  conjugate  gradient  method.  We  investigate 
different  pre-conditioners  in  the  nonlinear  conjugate  gradient  method  and  compare  their  ef¬ 
ficiency  with  the  commonly  used  Gauss-Newton  method.  By  using  a  pre-conditioner  based 
on  the  Hessian  matrix  of  the  objective  function,  the  nonhnear  conjugate  gradient  method 
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requires  fewer  iterations  than  the  Gauss-Newton  algorithna,  however  it  may  not  be  as  effi¬ 
cient  in  terms  of  CPU  time.  Simpler  pre-conditioners  require  more  iterations  but  less  CPU 
time  since  the  computation  of  the  Hessian  is  avoided. 

We  demonstrate  our  new  algorithms  on  a  large  dipole-dipole  data  set  collected  at  a 
groundwater  contamination  site  in  the  Aberjona  Watershed,  Woburn,  Massachusetts.  The 
objective  is  to  obtain  a  3-D  resistivity  model  of  the  subsurface  to  constrain  hydrological  mod¬ 
els  of  the  area.  With  this  data  set,  the  nonlinear  conjugate  gradient  algorithm  finds  minimum 
structure  models  with  approximately  ten  times  less  computation  than  the  Gauss-Newton 
method.  The  results  correlate  well  with  other  geophysical  results  at  the  site,  including  GPR 
sections  and  cone  penetrometer  logs. 

1.  Introduction 

The  d.c.  electrical  resistivity  method  is  a  principal  geophysical  exploration  technique  which 
has  been  used  extensively  for  subsurface  characterization.  The  implementation  of  the  method 
consists  of  injecting  current  into  the  ground,  measuring  the  electrical  potential  at  intervals  on 
the  surface  or  in  boreholes,  and  from  those  measurements  deducing  the  subsurface  resistivity 
distribution.  Because  the  instrumentation  is  simple  and  the  data  acquisition  is  straightfor¬ 
ward,  the  method  is  extremely  cost  effective.  Important  applications  include  mineral  deposits 
prospecting  (Keller,  1966j  Burger,  1992),  groundwater  exploration  (McNeill,  1990;  Medeiros 
and  Lima,  1990),  geothermal  studies  (Burger,  1992)  and  engineering  construction  projects 
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(Bogoslovsky  1979;  Shima,  1992).  Recently,  this  method  is  being  increasingly  applied  to 
environmental  protection  surveys  owing  to  its  non-invasive  property.  The  applications  in¬ 
clude  general  hydro-geological  mapping  {Okko,1993),  monitoring  contaminated  fluid  migra¬ 
tion  (Blum,  1989;  Van  et  al,  1992;  Spies  et  al,  1995),  mapping  the  extent  of  landflll  areas 
(Carpenter  et  al,  1990)  and  detection  of  contaminant  plumes  (Mazac  et  al,  1990a;  Buselli 
et  al,  1991). 

The  electrical  resistivity  of  geological  materials  depends  on  mineralogy,  clay  content, 
pore  fluid,  permeability,  conducting  metal  content,  and  other  properties  of  the  materials. 
The  resistivity  value  of  different  geological  materials  can  vary  from  10~^Q.m  (Pyrrhotite) 
to  10^'^n.m  (Dry  Limestone),  a  range  in  values  that  may  be  the  widest  of  any  common 
physical  property  of  earth  materials.  The  electrical  resistivity  method  is  ideally  suited  for  the 
task  of  detection  of  chemical  contamination.  The  resistivity  of  some  chemical  contaminants 
may  be  much  in  contrast  with  the  surrounding  natural  materials,  this  contrast  makes  the 
resistivity  method  the  most  sensitive  in  discerning  contaminants  from  groundwater.  Even 
when  the  resistivity  image  does  not  reveal  the  presence  of  chemical  contaminants  directly, 
it  can  provide  valuable  information  on  soil  properties  that  control  the  transport  of  chemical 
contaminants. 

The  task  of  the  d.c.  electrical  resistivity  inversion  is  to  solve  for  the  subsurface  resistivity 
distribution  from  the  measurements  of  electrical  potential  on  the  surface  of  the  earth  or  in 
boreholes.  This  is  a  diflicult  problem  on  many  levels.  One  major  difficulty  encountered  is 
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that  the  d.c.  resistivity  inverse  problem  is  ill-posed.  The  ill-posedness  comes  from  the  fact 
that  the  number  of  measurements  is  always  finite  while  the  unknown  resistivity  is  a  contin¬ 
uous  function  which  contains  in  principle  an  infinite  number  of  variables.  It  is  impossible 
to  construct  a  solution  that  is  stable  and  unique  based  on  fitting  data  alone.  Therefore 
solutions  must  incorporate  some  a  priori  information  or  be  regularized  in  order  to  be  physi¬ 
cally  meaningful  and  stable  against  the  noise  in  the  data.  Since  such  a  priori  information  or 
regularization  is  usually  difficult  to  obtain  directly  from  the  geological  reality,  it  is  subjected 
to  personal  bias.  Which  o  priori  assumption  or  regularization  criteria  is  thus  appropriate 
for  3-D  resistivity  inversion  remains  a  problem.  The  other  major  difficulty  associated  with 
the  inversion  is  a  result  of  the  problem’s  nonlinear  and  numerically  intensive  nature.  Since, 
in  practice,  a  moderate  3-D  resistivity  model  always  involves  a  large  number  of  model  pa¬ 
rameters  and  data  set,  the  forward  modeling  entails  solving  a  large  matrix  system.  Further, 
because  the  observed  electrical  potential  data  are  nonlinearly  dependent  on  the  resistivity 
parameters,  iterative  methods  are  needed  to  obtain  the  inversion  solutions.  Owing  to  this 
computational  difficulty,  methods  for  resistivity  inversion  were  mostly  restricted  to  1-D  (In¬ 
man,  1975;  Parker,  1984)  or  2-D  (Pelton  et  ai,  1978;  Tripp  et  ai,  1984;  Smith  and  Vozoff, 
1984;  Narayan,  1990)  earth  models  in  the  past  decades,  only  few  3-D  model  inversions  were 
found  in  recent  Uteratures  (Patrick  et  al,  1981;  Park  et  al,  1991;  Ellis  and  Oldenburg  1994; 
Zhang  et  al.,  1995). 

In  the  past,  most  methods  for  d.c.  resistivity  inversion  used  coarsely  parameterized  models 
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(e.g.,  large  layers  or  blocks)  to  make  the  problem  well-posed.  These  solutions  suppress  signif¬ 
icant  structures  and  can  hardly  match  the  geological  reality.  A  more  objective  approach  that 
allows  finely  parameterized  model  was  developed  by  Parker  (1984)  using  bilayer  expansion 
of  Green’s  function  on  a  1-D  uniform  grided  earth  model.  He  introduced  some  smoothness 
constrains  to  eliminate  the  rapid  unstable  oscillation  of  the  resistivity  values  that  are  caused 
by  the  non-uniqueness  of  the  problem.  Since  the  bilayer  expansion  was  formulated  based 
on  a  simple  layered  model,  his  approach  was  limited  to  1-D  earth  models.  Park  and  Van 
(1991)  and  Zhang  et  al.  (1995)  developed  a  3-D  resistivity  inversion  based  on  the  maximum 
likelihood  method  (Tarantola  and  Valette,  1982;  Madden,  1990)  using  statistical  information 
to  constrain  the  models.  The  method  philosophically  requires  full  statistical  knowledge  of 
the  model  parameters,  however,  in  practice,  such  information  is  often  unattainable,  there¬ 
fore  the  statistical  inference  was  replaced  by  a  uniform  damping  weighted  on  each  model 
parameters  which  makes  the  method  effectively  the  same  as  the  conventionally  damped  least 
square  method. 

The  approach  presented  in  this  report  is  based  on  Tikhonov  regularization  (Tikhonov, 
1962).  The  method  finds  solutions  by  emphasizing  the  importance  of  the  spatial  correlation 
of  the  model  parameters,  it  seeks  spatially  smooth  (or  “minimum  structure”  solutions  of 
the  inverse  problem.  The  approach  has  been  taken  by  Pilkington  and  Todoeschuck  (1992) 
in  cross-hole  seismic  tomography,  by  Scale  et  al.  (1990)  in  refraction  seismic  profile,  and  by 
Jiracek  et  al.  (1987)  and  deGroot-Hedlin  et  al.  (1990)  in  2-D  magnetotellurics.  These  studies 
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indicate  that  the  minimum  structure  approach  is  very  useful  in  the  case  of  highly  non-unique 
problems,  it  suppresses  model  structures  that  are  not  essential  in  matching  the  observation 
and  provides  a  lower  bound  on  model  complexity.  In  the  case  of  the  electrical  resistivity 
method,  we  believe  the  minimum  structure  models  are  reasonable  representations  of  the  real 
earth  because  the  diffusive  nature  of  the  electrical  energy  propagation  can  not  resolve  sharp 
resistivity  contrasts. 

When  miuiTTnim  structure  models  are  sought,  the  inversion  requires  insertion  of  some 
smoothness  constraints  which  are  used  to  minimize  the  model  roughness.  One  objective  of 
this  report  is  to  investigate  the  effect  of  the  smoothness  constraints  used  in  the  Tikhonov 
regularization  method.  We  will  compare  three  commonly  used  smoothness  constraints,  i.e., 
the  zeroth,  the  first  and  the  second  order  spatial  derivatives  of  model  parameters,  and  will 
show  that  their  effect  on  the  stability  of  the  solution  rely  on  the  behavior  of  the  sensitivity 
function,  which  is  a  measure  of  how  the  electrical  potential  outputs  change  due  to  a  small 
perturbation  in  the  resistivity  model  parameters. 

The  difficulty  of  implementing  such  minimum  structure  regularization  on  a  3-D  resistivity 
inversion  depends,  in  part,  on  the  numerical  algorithm  used  for  minimization.  Because  of 
the  nonlinearity  of  the  forward  problem,  iterative  minimization  schemes  are  needed  to  obtain 
solutions.  The  nonlinearity  may  also  induce  multiple  local  minima  in  the  objective  function 
used  to  find  mininumi  structure  models  that  fit  the  data,  it  is  therefore  necessary  to  repeat 
the  iterative  procedure  by  varying  initial  models  for  the  minimization  algorithm.  Further, 
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due  to  the  non-uniqueness  and  uncertainty  of  the  problem,  it  is  desirable  to  find  the  full  set 
of  acceptable  solutions,  or  at  least,  find  as  many  solutions  that  fit  the  data  as  possible.  In 
doing  so,  one  needs  to  repeat  the  iteration  algorithm  by  varying  the  a  priori  model  or  the 
smoothness  constraint.  Current  resistivity  inversion  algorithms  are  solved  by  an  iterative 
linearized  procedure  that  is  often  called  the  Gauss-Newton  method  (Tarantola  and  Valette, 
1982;  Park  and  Van,  1991;  Zhang  et  a/.,  1995).  It  starts  from  an  initial  model,  then  estimates 
a  small  perturbation  of  the  current  model  based  upon  the  Taylor  expansion,  this  process  is 
then  repeated  until  the  solution  converges.  This  method  suffers  from  slow  convergence  at 
the  early  stage  of  iterations.  It  has  been  shown  that  it  may  take  tens  of  hours  of  CPU  time 
on  a  high  speed  workstation  to  find  a  single  inversion  solution  on  a  small  (20x20x10)  model 
(EUis,  1995).  A  faster  and  robust  algorithm  is  very  desirable  for  extensively  studying  the 
uncertainty  and  the  resolution  of  the  problem  and  to  benefit  the  3-D  resistivity  method  in 
future  environmental  and  engineering  application. 

The  second  objective  of  this  report  is  to  investigate  a  more  efficient  algorithm  based 
on  the  conjugate  gradient  method  to  solve  the  nonlinear  minimization.  We  will  show  that 
with  an  appropriate  pre-conditioner  our  algorithm  requires  less  computer  resources  than  the 
Gauss-Newton  method  and  its  superiority  will  become  important  when  the  initial  model  is 
far  away  from  the  true  model. 

We  finally  demonstrate  the  inversion  algorithm  by  applying  it  to  actual  field  data  gathered 
using  a  dipole-dipole  electrical  resistivity  survey  at  a  superfund  contamination  site  in  the 
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Aberjona  Watershed,  Woburn,  Massachusetts.  The  primary  goal  of  that  experiment  is  to 
construct  a  3-D  subsurface  resistivity  structure  which  can  be  used  to  augment  hydrological 
models  of  the  area  in  an  effort  to  understand  the  transport  mechanisms  that  have  been 
involved  in  the  contamination  of  groundwater.  The  nonlinear  conjugate  gradient  algorithm, 
will  be  used  to  find  miniTmim  structure  models,  and  its  efficiency  will  be  compared  to  the 
Gauss-Newton  method.  Our  inversion  results  are  correlated  with  other  geophysical  results  at 
the  site  to  examine  their  consistency  and  accuracy,  and  are  found  to  be  in  good  agreement. 

2.  Theory 

2.1.  3-D  Electrical  Resistivity  Forward  Modeling 

The  d.c.  electrical  resistivity  forward  problem  is  the  solution  for  electrical  potential  in  an 
inhomogeneous  medium  governed  by  the  equation 

^  •  l~T~ - Vy  ^)]  = 

where  p{x,  y,  z)  is  the  resistivity  [fi.m]  of  the  medium,  j{x,  y,  z)  is  the  current  source  density 
[Am“^],  and  v{x,  y,  z)  is  the  electric  potential  [v]  subject  to  appropriate  boundary  conditions. 
On  the  surface  of  the  earth,  it  is  necessary  to  use  Neumann  boundary  condition,  ||  =  0, 
where  n  is  the  direction  normal  to  the  boundary.  On  portions  of  the  boundary  inside  the 
earth,  an  exact  boundary  condition  is  not  available  but  various  approximate  boundary  con¬ 
ditions  including  Dirichlet  and  mixed  boundary  conditions  (Day  and  Morrison,  1979;  Zhang 
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et  al,  1995)  can  be  used.  For  numerical  simplicity,  we  assume  that  the  model  boundaries 
are  far  from  the  source  and  receiver  so  that  a  Dirichlet  boundary  condition,  u  =  0,  can  be 
used. 

To  solve  (1)  numerically,  we  use  the  transmission  network  analog  developed  by  Madden 
(1972)  to  discretize  the  3-D  model  into  a  network  that  consists  of  network  nodes,  boundary 
nodes,  and  impedance  branches  (Figure  1).  Equation  (1)  is  then  approximated  by  a  linear 
system  of  equations: 

Kv  =  s  (2) 

where  v  is  a  vector  of  the  potentials  at  the  network  nodes,  s  is  the  current  source  vector,  and 
K  is  a  real,  symmetric,  and  positive-definite  matrix  which  depends  on  the  resistivities  and 
dimensions  of  the  network  cells.  Thus,  the  resistivity  function  p(x,  y,  z)  is  represented  by 
discrete  values  To  efficiently  solve  the  forward  problem  (2),  a  linear  conjugate  gradient 
algorithm  is  used  in  this  paper.  Details  of  this  solution  can  be  found  in  Zhang  et  al.  (1995). 

2.2.  Inversion  Method 

A  practical  d.c.  electrical  resistivity  inverse  problem  may  be  defined  as:  Given  a  set  of 
electrical  potential  measurements  d  =  (di,d2,  made  at  the  surface  or  in  boreholes, 

determine  as  much  information  as  possible  about  the  subsurface  resistivity.  This  may  be 
written  via  an  equation, 

d  =  Gtti  e  (3) 
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where  m  =  m[x)  is  the  unknown  resistivity  functional,  G  represents  the  forward  modeling 
operator  which  maps  the  model  space  to  the  data  space,  and  e  represents  the  error.  Since 
the  potential  measurements  can  be  made  only  at  a  finite  number  of  locations,  whereas  the 
resistivity  functional  contains  in  principal  an  infinite  number  of  variables,  this  problem  is 
ill-posed.  Solutions  constrained  only  by  the  data  alone  can  never  be  unique.  It  is  necessary 
to  incorporate  a  priori  information  in  order  to  define  a  unique  solution. 

Classical  remedies  that  advocate  a  priori  preference  into  the  model  fall  under  two  classes. 
One  class  assumes  that  the  model  parameters  are  random  variables  so  that  statistical  in¬ 
formation  can  be  introduced  to  constrain  the  model.  Important  approaches  include  the  fol¬ 
lowing:  the  Bayesian  inference  (Duijndam,  1988);  the  stochastic  inversion  (Franklin,  1970), 
and  the  maximum  likelihood  method  (Tarantola  and  Valette,  1982;  Madden,  1990).  The 
second  class  assumes  some  “regularity”  properties  of  the  solution  such  as  a  constraint  on 
the  spatial  smoothness  of  the  model  parameter.  This  technique  is  known  as  the  Tikhonov 
regularization  (Tikhonov,  1962).  It  is  important  to  note  that  the  methods  in  the  first  class 
philosophically  require  full  statistical  knowledge  of  the  model  parameters.  However,  often  in 
practice  when  such  information  is  unavailable,  these  methods  may  impose  some  smoothness 
constrains  which  makes  them  effectively  the  same  as  the  Tikhonov  regularization. 

Tikhonov  regularization  defines  a  solution  of  the  inverse  problem  that  fits  the  data  but 
also  that  has  minimum  possible  structures.  In  case  of  the  highly  non-unique  problems,  this 
technique  is  very  useful  because  it  provides  the  simplest,  or  “minimum  structure”  solutions 
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(Constable  et  al,  1987).  Solving  for  the  smoothest  model  can  avoid  being  misled  by  features 
that  appear  in  the  model  but  are  not  essential  in  matching  the  observation.  Therefore, 
anomalies  that  appear  in  the  model  can  only  come  from  the  data. 

Using  a  least-square  criterion,  Tikhonov  regularization  defines  a  solution  that  is  a  joint 
minimization  of  data  misfit  and  a  “stabilizing  functional”: 

$  =  (d  -  (d  —  Gm)  -f  r(m  -  mo)^L^L(m  -  mo)  =  min  (4) 

where  ^  is  the  objective  functional  to  be  minimized,  Rdd  is  data  covariance  matrix,  L  is  a 
linear  operator,  r  is  a  positive  number  known  as  the  regularization  parameter,  and  mo  is  a 
priori  model. 

The  first  term  in  the  objective  function  measures  the  data  misfit.  It  is  the  likeUhood 
function  when  noise  contaminating  each  observed  data,  e,  has  Gaussian  normal  distribution 
with  zero  mean  and  variance  The  second  term  defines  the  stabilizing  functional  as  the 
regularization  term  which  measures  the  spatial  roughness  of  squared  norm  of  L(m  -  mo). 
In  stochastic  or  the  maximum  hkelihood  inversion,  m©  is  taken  to  be  a  priori  mean  of  m 
and  L  is  chosen  such  that  (L^L)“^  =  Rmm  is  a  a  priori  covariance  of  m.  In  a  miuinrmm 
structure  approach,  L  is  a  diflEerential  operator  and  mo  is  taken  to  be  a  simple  priori  model. 
In  practice  it  is  desirable  to  vary  mo  so  that  multiple  solutions  can  be  obtained. 

Conventionally,  the  stabilizing  functional  has  three  options: 

1,  L  = /,  (m  —  mo)^L^L(m  —  mo)  =  j dxdydz\m—mQ\^  (5) 
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(6) 


2,  L  =  V,  (m  -  mo)^L^L(m  -  mo)  =  j  dxdydz\V {m  -  mo)P 

3,  L  =  V^  {m-mQ)VL{m-mo)  =  J dxdydz\V^{m-mo)\‘^.  (7) 

It  is  not  certain  whether  any  of  those  options  guarantees  a  well-posed  minimization  in  3-D 
resistivity  inversion.  In  the  next  section  we  will  discuss  this  issue. 

3.  Comparison  of  Stabilizing  Functionals 

Among  proponents  of  the  minimum-structure  approach,  there  is  no  consensus  on  the  best 
smoothing  operator  L  to  use.  Constable  (1990)  defined  L  in  terms  of  the  first  derivative, 
as  in  Eq  (6).  Scales  et  al.  (1990)  defined  L  in  terms  of  the  second  derivative,  as  in  Eq  (7). 
Ellis  and  Oldenburg  (1994)  defined  L  as  a  combination  of  both  the  first  and  the  second 
derivatives.  To  investigate  their  influences  on  the  3-D  resistivity  problem,  we  compare  them 
theoretically  and  numerically. 

We  know  that  in  Eq  (4),  when  ^  is  minimized  its  first  order  partial  derivative  with 
respect  to  m  is  zero  yielding, 

(d  -  Gm)  +  rL^L(m  -  mo)  =  0  (8) 

where  m  is  the  solution  at  which  ^'(m)  is  minimized.  A  is  the  sensitivity  operator,  or  Frechet 
derivative.  For  a  discrete  model,  A  is  the  matrix 


which  measures  how  data  changes  at  the  ith  receiver  to  a  change  of  jth  resistivity.  Park  et 
al.  (1991)  derived  that  the  sensitivity  matrix  is  given  by  the  inner  product  of  the  current 
density  ( J^)  from  a  point  source  at  the  transmitter  and  the  current  density  {Jr)  from  a  point 
source  at  the  receiver  integrated  over  the  perturbed  volume, 

Ay  j  Js'  (10) 

In  the  3-D  problem,  since  the  current  distribution  from  a  point  source  approaches  infinity  at 
the  source,  this  sensitivity  function  has  singularities  at  the  location  of  source  and  receivers. 
Any  small  perturbation  in  the  data  will  result  in  a  large  variation  in  the  model  space.  This 
aspect  of  the  physics  must  be  accounted  for  in  the  regularization  method. 

Using  Eq  (8),  we  obtain  the  following  relationship  between  the  regularized  solution  model 
m  and  the  sensitivity  matrix  A, 

L^Lm  =  L^Lmo  +  r-i^AoRji(d-Gm)  (11) 

i 

Therefore,  when  L  is  an  identity,  the  solution  model  m  equals  mo  plus  a  linear  combination 
of  the  sensitivity  matrix  multipUed  by  the  data  residue.  When  L  is  the  first  order  spatial 
derivative  operator,  the  Laplacian  of  7n  is  a  linear  combination  of  the  sensitivity  matrix 
multiplied  by  the  data  residue.  When  L  is  the  second  order  spatial  "derivative  operator,  the 
Laplacian  squares  of  m  is  a  linear  combination  of  the  sensitivity  matrix  multiplied  by  the 
data  residue.  To  test  which  choice  of  the  stabihzing  functional  yields  a  stable  solution,  we 
design  a  simple  synthetic  test  problem. 
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The  model  is  a  conductive  block  (in.m)  buried  in  a  homogeneous  half  space  (lOOfi.m). 
The  model  is  discretized  into  21x21x15  elements  with  10  m  spacing.  The  (Ifi.m)  conductive 
block  is  discretized  into  7x7x3  elements.  25  receiver  electrodes  are  placed  on  the  surface, 
among  which  nine  of  them  are  also  used  as  current  electrodes.  A  total  of  (9x24=216) 
observations  are  produced  by  forward  modeling  and  3%  random  Gaussian  noise  is  added  to 
the  data. 

Before  investigating  the  three  stabihzing  functionals,  the  amount  of  smoothing  done 
by  each  stabilizing  functional  needs  to  be  justified.  Previous  studies  (Backus  and  Gilbert, 
1970;  Parker,  1984)  indicate  that  when  noise  contaminates  the  data  there  is  a  trade-off 
between  the  data  misfit  and  the  model  roughness.  By  varying  the  regularization  parameter 
r,  a  trade-off  curve  (Figure  3a)  is  generated  using  stabilizing  functional  (7)  as  an  example. 
Three  points — ri,  T2,  and  tz — are  selected  form  the  curve  and  their  associated  models  are 
reconstructed  and  shown  in  Figure  3b.  It  is  found  that  when  the  model  is  too  smooth  (when 
r  =  5.5),  it  cannot  fit  the  data  very  well  (x^  =  275).  As  r  is  lowered  (r  =  1.1),  the  model 
has  a  smooth  structure  and  provides  a  better  reconstruction  (x^  =  176).  Over-fitting  data 
(x  =  130  at  r  =  0.1)  results  in  a  rough  model  which  contains  incorrect  surface  anomalies. 
In  general,  if  data  contains  the  random  Gaussian  noise,  the  expected  value  of  x^  equals  to 
the  number  of  independent  data,  therefore  choosing  to  fit  to  the  number  of  independent 
data  is  an  optimum  choice. 

In  our  case,  we  fit  the  potential  data  to  x^  misfit  of  180  (number  of  independent  source  and 
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receiver  pairs).  The  results  for  an  inversion  obtained  using  the  zeroth,  the  first  and  the  second 
order  regularization  are  shown  in  Figure  4.  The  one  with  the  zeroth  order  regularization 
has  large  resistivity  variations  in  the  vicinity  of  source  and  receivers  indicating  that  the 
singularities  in  the  sensitivity  functions  are  not  suppressed  by  this  stabilizing  functional. 
The  first  order  regularization  yields  better  results,  however  the  surface  artifacts  are  still  seen 
in  the  resulting  model.  The  second  order  regularization  successfully  suppresses  the  surface 
artifacts,  giving  the  smoothest  result. 

4.  Comparison  of  Minimization  Algorithni 

In  the  preceding  sections  we  have  defined  the  solution  of  the  resistivity  inverse  problem  by 
minimizing  the  objective  function  'F  in  equation  (4).  Since  the  forward  modeling  operator 
G  depends  nonlinearly  on  m,  'F  is  non-quadratic  and  an  iterative  minimization  is  required. 

In  this  section  we  will  investigate  two  algorithms  for  minimizing  the  objective  functional, 
i.e.,  the  Gauss-Newton  method,  and  a  nonlinear  conjugate  gradient  method. 

4.1.  Gauss-Newton  Method 

The  Gauss-Newton  method  is  based  on  expanding  G  in  a  Taylor  series  and  calculating  the 
model  correction  at  each  iteration  on  the  assumption  of  local  linearity.  By  Taylor  expansion, 
we  have 

G(m  +  6m)  =  G(m)  4-  A6m  (12) 
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where  we  have  ignored  higher  order  derivatives.  Thus  the  value  of  ^  predicted  by  (12)  is 

^(m  +  6m)  =  (d  -  Gm  -  A6Tn)^R7j(d  -  Gm  —  A6m)  + 

r(m  -  mo  +  6m)^L^L(m  -  mo  +  6m)  (13) 

With  this  approximation,  depends  quadratically  on  6m  and  is  minimized  by  setting 
d^/d6m  =  0,  thus  6m  is  found  by  solving 

6m  =  (A^R^^A  +  rL^L)-^(A^RJdH<^  -  Gm)  +  rL^L(mo  -  m))  (14) 

The  Gauss-Newton  method  thus  constructs  a  sequence  of  the  models  by 

mfc+i  =  m/fc  4-  6mk  (15) 

where  6mk  solves  the  linear  equation 

6mk  =  (Ak^R^^Ak  +  rL^L)-HAk^Rji(d  -  Gm^)  +  TL^L(mo  -  m*))  (16) 

and  a  repetition  of  this  process  yields  successive  estimates  mi,  m2,  ...  mjt  until  the  minimum 
is  found. 

Solving  this  system  by  computing  the  inverse  of  the  Hessian,  (A^R^J  A-l-rL^L),  requires 
a  tremendous  amount  of  computing  resources.  Zhang  et  al.  (1995)- suggested  that  one  can 
solve  for  6m  without  direct  computation  of  the  Hessian  matrix  by  applying  a  linear  conjugate 
gradient  scheme.  This  approach  significantly  reduces  the  amount  of  computing  resources  and 
makes  the  3-D  resistivity  inversion  using  the  Gauss-Newton  method  practical. 
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The  Gauss-Newton  method  has  the  desirable  feature  of  rapid  convergence  if  the  initial 
guess  is  close  to  the  solution,  especially  if,  in  the  neighborhood  of  the  solution  the  objective 
function  appears  to  become  quadratic,  then  Gauss-Newton  method  will  find  the  solution  in 
one  step.  Unavoidably,  an  initial  guess  that  is  far  from  the  solution  may  lead  to  difficulty. 
The  convergence  may  reveal  overshooting  and  unstable  features.  Therefore  in  practice  a 
small  damping  is  often  added  to  the  Hessian  in  order  to  secure  the  stability.  The  penalty  is 
the  slow  convergence  at  the  early  iterations. 

4.2.  Nonlinezir  Conjugate  Gradient  Method 

Another  potentially  more  efficient  algorithm  for  minimizing  '^{m)  is  developed  in  this  paper. 
This  algorithm,  the  Nonlinear  Conjugate  Gradient  (NLCG),  was  first  extended  by  Fletcher 
and  Reeves  (1959)  from  the  linear  conjugate  gradient.  This  method  yields  an  alternative 
minimization  of  a  nonlinear  function  directly  without  making  any  assumptions  about  its 
linearity. 

The  method  of  NLCG  is  outUned  as  follows: 

(a)  Choose  mi,  and  set  gi  =  -V'5'(mi),  hi  =  Cgi 

(b)  Find  a  =  cxk  which  minimizes  +  <^hk) 

(c)  nik+i  =mk  +  akhk,  gk+i  =  -V^(mfc+i) 

(d)  hk+i  =  Cgk+i  -I-  Pkhk,  where  pk  = 

In  the  algorithm,  gk  is  the  gradient  of  the  objective  functional,  hk  is  the  direction  along  which 
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the  parameter  ak  is  searched  to  minimize  "^{nik-hoihk),  C*  is  a  preconditioning  operator  which 
we  will  discuss  later,  and  the  formula  for  hk  and  gk  above  imply  the  basic  CG-relation 

Qi  ■  hj  =  0,  j)  (17) 

There  are  two  important  issues  of  this  algorithm  that  need  to  be  addressed.  First,  in  step 
(b)  there  is  a  line  minimization  required  so  as  to  find  a  single  variable  a  which  minimizes 
^(rrik  +  Oihk).  At  the  earlier  stage  of  the  iterations,  the  surface  of  ^  tends  to  be  more 
distorted  due  to  a  relatively  large  amount  of  nonlinearity  of  G,  while  at  the  later  stage  of 
the  iterations  it  becomes  more  quadratic  in  the  vicinity  of  the  minimum.  This  behavior 
of  'i'  leads  us  to  design  a  line  minimization  routine  that  is  adequate  for  both  the  earlier 
and  the  later  iterations.  The  strategy  is:  at  every  iteration  k,  try  the  first  step  of  the  line 
minimization  based  on  a  quadratic  assumption,  a  =  After  a  minimum  is  bracket 

then  estimate  the  minimum  point  and  estimate  its  value  by  using  a  cubic  interpolation.  Such 
a  strategy  would  not  only  work  at  the  earlier  stage  of  iteration  but  also  guarantee  a  rapid 
quadratic  convergence  at  the  later  stage  of  the  iterations. 

Second,  because  poor  conditioning  of  ^  leads  to  slow  convergence,  it  is  necessary  to  incor¬ 
porate  a  preconditioner  into  the  algorithm.  A  good  preconditioner  may  help  to  decrease  the 
objective  function  more  quickly  in  the  early  iterations.  If  the  problem  were  linear  (quadratic 
),  a  perfect  preconditioner  would  be  the  inverse  of  the  Hessian  matrix 

C  =  A  -I-  rL^L)-'  (18) 
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It  is  probably  safe  to  assume  that  this  preconditioner  works  well  in  the  nonlinear  problem, 
too.  However,  applying  it  requires  extensive  computing  resources.  Hence,  we  explore  an 
improved  preconditioner  that  approximates  (18)  and  is  fast  to  implement.  We  suggest  the 
use  of  the  following  formula  as  our  preconditioner, 

C  =  {tL^L  +  sI)-^  (19) 

where 

(20) 

v/u 

in  which 

u  =  [logpo,  logpo,  logpo,  ....logpo]  (21) 

is  a  model  with  constant  resistivity  po- 

The  substitution  of  (19)  for  (18)  scales  the  d.c.  component  of  the  term  A^RjjA.  Nu¬ 
merically,  the  constant  s  defined  by  Eq  (20)  is  easy  to  calculate,  because  according  to  Eq  (9), 
we  have, 

(Au)i  =  {Gm)i  (22) 

which  is  just  the  response  of  the  constant  model  and  can  be  obtained  with  no  computation 
beyond  the  forward  calculation. 

One  important  issue  needs  to  be  addressed.  Because  of  the  strong  nonlinearity  of  the 
resistivity  problem,  the  objective  function  ^  may  has  not  only  a  global  minimum  but  also 
some  undesired  local  minimums.  Although  both  the  Gauss-Newton  and  NLCG  methods 
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solve  the  nonlinear  minimization  problem  and  will  find  a  minimum,  they  can  not  guarantee 
if  the  TniniTmiTn  that  is  found  is  the  global  minimum.  Finding  the  global  minimum  is,  in 
general,  a  very  difficult  problem.  A  widely  used  method  is  to  repeat  the  algorithm  starting 
from  different  initial  models,  finding  minimums,  and  then  picking  the  most  extreme  of  them 
(if  they  are  not  all  the  same)  (Press,  1992).  A  faster  algorithm  will  be  benefit  for  this 
procedure. 

4.3.  Numerical  Comparison 

We  conducted  a  number  of  convergence  tests  for  the  Gauss-Newton  and  the  nonlinear  con¬ 
jugate  gradient  method  using  the  same  synthetic  model  described  before  (Figure  2).  The 
Gauss-Newton  minimization  scheme  used  in  this  report  was  developed  by  Zhang  6t  al.  (1995). 

Choosing  a  starting  model  with  a  homogeneous  resistivity  value  of  200  Q.m,  the  conver¬ 
gence  results  for  each  optimization  method  are  shown  in  Figure  5a.  Here  ‘NLCG’  denotes 
the  NLCG  method  with  the  improved  preconditioner  (19),  ‘NLCG-1’  denotes  the  NLCG 
method  without  preconditioning  (C=identity),  and  ‘NLCG-2’  denotes  the  NLCG  with  the 
Hessian  (18).  From  these  results  it  is  clear  that  the  convergence  rate  of  the  NLCG  method 
depends  strongly  on  the  preconditioner:  NLCG  is  the  most  efi&cient  method  in  terms  of  the 
CPU  time,  and  NLCG-2  is  the  least  efficient  method.  We  also  see  that  NLCG-2  makes  larger 
decreases  in  ^  but  takes  more  CPU  time.  The  Gauss-Newton  method  interpolates  between 
NLCG  and  NLCG-2.  When  all  of  the  methods  converge,  they  find  the  same  solution  (Figure 
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5b). 


If  the  initial  model  is  chosen  to  be  close  to  the  true  model:  mo  :=  SOfi.m  {inside  the 
conductive  block);  =  lOOfi.m  {other  region),  the  Gauss-Newton  method  will  be  as  efficient 
as  that  of  the  NLCG  method  (Figure  6a).  This  result  demonstrates  that  the  Gauss-Newton 
method  has  a  fast  convergence  rate  in  the  neighborhood  of  the  solution.  Both  methods  find 
the  same  solution  (Figure  6b). 

When  the  initial  model  is  far  from  the  solution,  mo  =  50012. m,  the  NLCG  method  has 
much  faster  convergence  than  the  Gauss-Newton  method  (Figure  7a).  Both  methods  finally 
find  the  same  solution  (Figure  7b). 

5.  Application  to  Field  Data:  3-D  Resistivity  Survey  at  the  Aber- 
jona  Watershed,  Woburn,  Massachusetts 
Background  and  Objective 

The  Aberjona  Watershed  is  located  in  eastern  Massachusetts  (Figure  8,  McBrearty,  1993). 
The  area  has  a  history  of  industrial  contamination  dating  back  to  the  beginning  of  the  cen¬ 
tury.  Industries,  including  leather  tanners,  metal  cleaning  and  automobile  salvage  yards, 
dumped  their  waste  into  the  ground  and  contaminated  the  ground  water  with  high  concen¬ 
trations  of  chemicals  shown  to  be  cancer-causing  in  laboratory  animals,  including  chloroform, 
perchloroethylene  (PCE),  trichloroethane,  and  trichloroethylene  (TCE).  A  suspected  result 
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of  the  contamination  was  a  high  incidence  of  childhood  leukemia  in  the  area  during  the  mid 
1980s  (Diperna,  1985).  The  area  was  designated  as  a  U.S.  Environmental  Protection  Agency 
superfund  site  in  1982,  and  has  become  the  focus  of  numerous  environmental  investigations. 

The  subsurface  geology  of  the  Aberjona  Watershed  is  characterized  by  glacial  outwash  and 
till  deposits.  One  of  the  major  valleys— the  Aberjona  valley  formed  by  the  glacial  movement, 
runs  beneath  the  main  branch  of  the  Aberjona  river.  The  valley  that  forms  the  transmissive 
regions  of  the  aquifer  to  a  depth  of  over  100  feet  is  filled  with  glacial  outwash  consisting  of 
fine  to  medium  sands  with  traces  of  silt  at  shallow  depth,  and  medium  to  coarse  sand  and 
gravel  as  the  depth  increases.  Previous  studies  also  indicated  a  peat  layer  overlaying  the 
stratified  sand  deposits  extending  over  a  large  area  with  the  thickness  range  from  2  to  7  feet 
(USGS  Report,  1989).  Because  of  the  shallow  water  table  the  soils  are  generally  saturated 
to  the  surface,  the  effect  of  the  topography  is  to  cause  a  damp  soft  soil  with  scattered  areas 
of  very  shallow  pooling  water.  Further,  during  the  winter  time,  the  hydraulic  conductivity 
decreases  causing  a  larger  areas  of  flooding. 

Our  study  site  (marked  inside  a  circle  in  Figure  8)  is  located  near  the  Well  H  region  of  the 
Aberjona  River.  Our  interest  focuses  on  characterizing  the  soil  around  the  river  in  an  effort 
to  understand  the  transport  mechanisms  that  have  been  involved  in  the  contamination  of  the 
drinking  water.  Before  our  resistivity  survey,  other  geophysical  investigations  were  employed 
in  the  area  including  ground  penetrating  radar,  and  cone  penetrometer  surveys  (Zeeb  et  a/., 
1994).  The  ultimate  goal  of  the  resistivity  survey  was  to  correlate  soil  properties  with  those 
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investigations  to  extend  the  point-source  ground  truth  information  into  a  3-D  map  of  the 
stratigraphy  of  the  whole  region. 

5.1.  The  Resistivity  Measurement 

We  performed  the  electrical  resistivity  survey  at  the  site  in  March  1995.  Figure  9  shows  the 
experiment  site  and  the  electrode  configuration.  There  were  80  electrodes  placed  at  25  ft 
intervals  in  an  8x10  array  on  the  surface,  covering  a  175  x  225  (/t^)  area.  Electrical  current 
was  injected  one  at  a  time  into  30  of  the  electrodes  (marked  in  red),  with  the  negative  end 
of  the  current  placed  1000  feet  south  of  the  Well  H.  Potential  differences  were  measured 
between  each  of  the  remaining  electrodes  and  a  point  adjacent  to  the  current  electrode.  A 
total  of  30x80  potential  differences  were  thus  obtained. 

Figure  10  shows  a  part  of  data  set  where  the  potential  differences  are  displayed  in  form  of 
contour  lines.  The  current  sources  are  clearly  noticeable  in  the  plot  where  they  correspond  to 
the  highest  potential  differences.  We  estimated  that  the  measured  potential  differences  were 
accurate  to  within  5  percent  based  upon  repeated  measurements,  such  that  the  standard 
deviation  for  each  datum  was  assumed  to  be  5  percent  of  the  potential  difference. 

5.2.  The  Inversion  Results 

We  generated  a  3-D  model  with  26x24x20  cells  to  cover  the  experimental  region.  To  construct 
a  starting  model,  surface  resistivity  information  was  used.  As  one  can  see  from  Figure  9,  a 
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part  of  the  Aberjona  River  meanders  its  way  southward,  about  100  feet  west  of  the  Well  H. 
The  river  was  built  into  the  starting  model  by  assigning  an  initial  resistivity  value  of  10  fi.m 
to  cells  traversed  the  path  of  the  river.  In  the  rest  of  the  area  resistivity  values  of  surface 
soil  samples  were  measured  and  their  mean  p  =  65(fi.m)  was  assigned  to  the  remaining  cells 

of  the  starting  model. 

The  NLCG  method  is  used  to  invert  the  3-D  resistivity  structure.  Gauss-Newton  is  also 
used  as  a  comparison.  The  convergence  rate  of  both  method  is  shown  in  Figure  11.  The 
Gauss-Newton  method  requires  approximately  15  hours  of  computing  time  on  a  DEC  alpha 
workstation  to  find  one  solution,  while  the  NLCG  methods  takes  1.5  hours. 

Figures  12a  and  12b  show  a  3-D  resistivity  inversion  model  as  horizontal  -I-  vertical 
cross-sections.  The  results  depict  a  resistive  zone  at  20  feet  below  the  surface,  dipping  from 
southeast  to  northwest.  Above  the  resistive  zone,  the  results  show  a  bowl-shape  conductive 
zone.  In  order  to  interpret  the  physical  soil  types,  a  vertical  cone  penetrometer  profile 
obtained  at  the  place  that  is  75  feet  west  of  the  Well  H  is  used  to  correlate  the  resistivity 
results.  The  results  are  shown  in.  Figure  13. 

The  cone  penetrometer  response  indicates  somewhat  complicated  soil  types  at  shallow 
depth.  Except  for  a  region  of  organic  peat  layers  directly  below  the  surface,  there  is  a  sec¬ 
tion  comprised  of  brown,  spongy  material  that  has  a  high  water  content  (80%  water  by 
weight)  and  contains  almost  no  organic  matter.  This  section  is  found  to  be  the  diatoma- 
ceous  earth  which  consists  primarily  of  siliceous  skeletons  of  diatom-microscopic.  Because 
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of  high  porosity  (70-90  %  in  general)  resulting  in  high  water  content,  the  diatomaceous 
earth  section  appears  to  be  very  conductive.  Adjacent  to  the  bowl-shaped  lower  boundary 
of  the  diatomaceous  earth  is,  as  expected,  a  glacial  sand  layer  which  is  quite  resistive.  The 
boundary  between  the  conductive  zone  (peats  -f  diatomaceous  earth)  and  the  resistive  zone 
(sand),  as  indicated  by  the  3-D  resistivity  inversion  model,  is  strongly  coherent  with  the 
cone  penetrometer  results.  However,  due  to  the  diffusive  nature  of  the  electrical  energy  and 
perhaps,  due  to  the  “minimum  structure”  approach,  the  3-D  resistivity  inversion  cannot 
delineate  the  heterogeneity  within  the  conductive  zone,  i.e.,  different  thin  peat  layers  and 
the  diatomaceous  earth  section. 

The  resistivity  results  are  in  strong  agreement  with  GPR  (Ground  Penetrating  Radar) 
experiments  (Cist  et  ai,  1995).  Figures  14a  and  14b  each  compare  a  vertical  GPR  section 
with  a  correspondent  vertical  cross-section  through  the  3-D  resistivity  inversion  model.  Prom 
these  plots,  one  can  see  that  the  shape  of  the  boundary  between  the  diatomaceous  earth  and 
the  sand  layer  is  in  good  agreement  between  the  two  experiments.  Below  this  interface,  the 
GPR  section  is  dominated  by  ringing  effects  of  the  antennas,  whereas  the  resistivity  model 
detects  vertical  and  lateral  variations  of  the  resistivity  distribution.  This  result  indicates  that 
the  resistivity  method  has  great  advantages  in  resolving  structures  at  depths  unattainable 
to  either  GPR  method  or  cone  penetrometer. 
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6.  Conclusions 


In  this  paper,  we  have  presented  a  method  for  solving  the  nonlinear  d.c.  electrical  resistivity 
inverse  problem.  To  deal  with  the  ill-posedness  of  the  problem,  the  method  uses  Tikhonov 
regularization  to  obtain  a  solution  which  has  the  least  structure  necessary  to  fit  the  data 
to  specified  error  bonds.  We  have  investigated  three  different  smoothness  constraints  and 
concluded  that  it  is  necessary  to  constrain  second  order  spatial  derivatives  of  the  resistivity 
function  to  obtain  smooth,  stable  solutions.  This  result  can  be  understood  in  terms  of  the 
sensitivity  functions  for  3-D  resistivity  data  and  the  fact  that  they  are  singular  at  source 
and  receiver  positions. 

To  implement  Tikhonov  regularization  numerically,  we  have  developed  a  fast  and  efficient 
algorithm  based  on  conjugate  gradient  with  preconditioning  to  solve  the  nonlinear  minimiza¬ 
tion  problem.  Using  a  simple  preconditioning  operator,  the  NLCG  algorithm  results  in  a 
tremendous  time  savings  over  the  conventional  Gauss-Newton  approach,  especially  when  the 
initial  guess  is  far  from  the  solution. 

We  have  successfully  applied  our  3-D  resistivity  inversion  algorithm  to  a  actual  resistivity 
field  survey  in  the  Aberjona  Watershed,  Woburn,  Massachusetts.  With  this  field  data,  the 
nonlinear  conjugate  gradient  algorithm  finds  minimum  structure  models  with  approximately 
ten  times  less  computation  than  the  Gauss-Newton  method.  Our  3-D  resistivity  image 
indicates  two  major  layers  of  the  site:  the  conductive  layer  which  involves  the  peat  and 
the  diatomaceous  earth  is  located  beneath  the  surface;  the  resistive  sand  layer  is  underlined 
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at  an  elevation  of  20  feet  deep,  dipping  from  south-east  to  north-west.  Results  have  been 
compared  and  found  to  be  strongly  correlated  with  other  geophysical  experiments  of  the  site 
such  as  cone  penetrometer  and  GPR  measurements. 
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0  -  network  node 

O  -  boundary  node 


Figure  1:  The  3-D  resistivity  model  is  discretized  into  a  network  which  consists  of  network 
nodes,  boundary  nodes,  and  impedance  branches. 
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View  of  the  Model 


current  receiver 


Figure  2;  The  synthetic  model  consists  of  the  7x7x3  prism  (p=l  ohm.m)  inside  a  21x21x15 
background  (p=100  ohm.m).  The  data  are  defined  on  a  5x5  electrode  array  and  comprise 
24  potential  field  measurements  for  each  of  9  pole  current  sources. 
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Figure  3a:  This  trade-off  curve  shows  how  the  data  misfit  is  balanced  by  the  model  roughness. 
If  the  model  is  too  smooth,  it  cannot  fit  the  data  well.  Over-fitting  the  data  will  result 
in  a  rough  model  which  contains  incorrect  surface  anomalies.  An  optimum  choice  is  to 
fit  chi-square  to  the  number  of  independent  data  which,  in  this  case,  is  180,  indicated  by 
the  dashed  line. 


Min.  Laplacian  Solutions  Controlled  by  x 
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Solutions  by  Three  Different  Stabilizing  Functionals 
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Efficiency  Test  of  Different  Methods 


Figure  5a:  NLCG  is  the  most  efficient  method  among  the  others.  NLCG-2  is  the  least 
efficient  method.  The  Gauss-Newton  method  interpolates  between  them.  The  initial 
model  for  this  test  is  200  (ohm.m). 

Solutions  by  Different  Algorithms  (m1=200  ohm.m) 
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Figure  6a:  When  the  initial  guess  is  close  to  the  true  model,  the  Guass-Newton  method  is 
about  as  efficient  as  the  NLCG  method.  Where  S  denotes  the  objective  function,  and 
x~2  denotes  the  chi-square.  The  initial  moel  for  this  test  is  ml=l-  (background),  =50 
(inside  the  conductive  block). 

Solutions  by  Different  Algorithms  (m1=  100(out),  50(in)) 


Gauss-Newton  NLCG 
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Figure  6b:  Gauss-Newton  finds  the  same  solution  as  that  of  NLCG  with  a  good  initial  model. 
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Figure  7a:  The  NLCG  method  is  much  more  efficient  in  terms  of  the  CPU  time  compared 
to  the  Gauss-Newton  method.  Where  S  denotes  the  objective  function,  X''2  denotes 
chi-square.  The  initial  model  is  cosen  as  m0=500  ohm-m. 

Solutions  by  Different  Algorithms  (m1=500) 
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Figure  7b:  Gauss-Newton  finds  the  similar  solution  as  that  of  NLCG  with  bad  intial  model. 
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BASE  MAP  MODIFIED  FROM  TOPOGRAPHICAL 
MAPST-12,  13.  16.  17.21.22.  CITY  OF  WOBURN. 
MASSACHUSETTS 


^  Wetland 
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Figure  8:  The  Aberjona  Watershed  is  located  in  eastern  Massachusetts.  Our  experiment  site 
(squared)  is  located  near  Well  H  region  of  te  Aberjona  river  in  the  town  of  Woburn,  a 
suburb  10  miles  north  of  Boston. 
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Figure  9:  The  resistivity  experiment  deployed  80  electrodes  deployed  in  the  Well  H  area,  30 
of  them  were  used  as  the  current  electrodes.  The  negative  current  electrode  is  1000  ft 
south  of  Well  H.  The  potential  differences  were  measured  with  respect  to  a  point  which 
is  adjacent  to  the  current  electrode. 
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Contour  of  the  observed  potential  differences 


-1.0  0.0 

Potential  Difference  (v) 


Figure  10:  Contour  plots  of  the  measured  potential  differences  from  four  current  sources. 
The  current  sources  are  noticeable,  which  appear  to  be  the  singularities. 
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Objective  Function  vs.  CPU  Time 


CPU  Time  (min.) 


Figure  11:  Iteration  histories  for  the  NLCG  and  Gauss-Newton  algorithm  of  the  Aberjona 
data.  The  NLCG  algorithm  yields  an  acceptable  inversion  model  approximately  ten  times 
faster  than  the  Gauss-Newton  method. 
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3-D  R9sislivity  Invereion  Model 
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Figure  12-(b)  3D  resistiuity  inuersion  model  (horizontal  slices)  for  the  Rberjona 
contamination  site.  The  resistiue  zone  at  depth,  dipping  from  east  to  west,  is 
interpreted  to  be  a  sand  formation  underlying  more  conductiue  diatomaceous  earth. 
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Figure  13,  Correlating  the  3-D  resistivity  inversion  model  with  the  Cone  penetrameter  resopnse. 
The  conductive  zone  is  found  to  be  a  region  of  peat  layers  and  a  diatomaceous  earth  section,  while 
the  resistive  zone  corresponds  to  a  glacial  sand  layer.  Note,  the  bowl-shape  lower  boundary  of  the 
diatomaceous  earth  is  in  good  agreement  between  the  two  methods. 
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Radar  Measurement  ( 75W) 


Resistivity  Measurement  ( 75W) 


<z> 


South  North 

Fig  14  -(a)  Comparison  of  resistivity  inversion  model  with  a  GPR  section 
located  at  75  ft  west  of  the  Well  H.  The  shape  of  the  boundary  between 
the  diatomaceous  earth  and  the  sand  layer  is  in  good  agreement. 
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Fig  14  -(b)  Comparison  of  resistivity  inversion  model  with  a  GPR  section 
located  at  150  feet  north  of  the  Well  H.  The  shape  of  the  boundary  between 
the  diatomaceous  earth  and  the  sand  layer  is  in  good  agreement. 
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ELECTROSEISMIC  WAVES  FROM  POINT  SOURCES  IN 

LAYERED  MEDIA 


Summary 

In  sedimentary,  porous  material  saturated  by  a  fluid  electrolyte,  mechanical  and  electromag¬ 
netic  disturbances  are  coupled.  The  coupling  is  electrokinetic  in  nature  due  to  an  excess 
of  electrolyte  ions  that  exists  in  a  fluid  layer  near  the  grain  surfaces  within  the  material. 
Seismic  waves  in  sedimentary  material  generate  relative  fluid-solid  motion  that  induces  an 
electrical  streaming  current.  When  a  seismic  pulse  traverses  contrasts  in  electrical  and/or 
mechanical  medium  properties  a  dynamic  streaming  current  imbalance  is  induced  across 
such  a  contrast,  generating  electromagnetic  disturbances  that  are  measurable  at  the  earth’s 
surface.  This  paper  numerically  determines  this  full-waveform  electroseismic  pointsource 
response  in  a  stratified  porous  medium.  It  is  shown  that  the  macroscopic  governing  equa¬ 
tions  controlling  the  coupled  electromagnetics  and  acoustics  of  porous  media  decouple  into 
two  sets  corresponding  to  vertical  a  (PSVTM)  or  horizontal  (SHTE)  polarization  of  the 
transverse  wavefields.  The  seismic  pulse  is  shown  to  induce  electric  fields  that  travel  with 
the  compressional  wavespeed  and  magnetic  fields  that  travel  with  the  shear  waves.  The 
frequency  content  of  the  converted  electromagnetic  field  has  the  same  frequency  content  of 
the  driving  incident  seismic  pulse,  as  long  as  the  propagation  distances  are  much  less  than 
the  electromagnetic  skin  depth. 
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Snapshots  in  time  and  converted  electromagnetic  amplitudes  versus  seismic  point  source- 
antenna  offset  are  calculated  for  contrasts  in  mechanical  and/or  electrical  medium  properties. 
The  TM  component  amplitude  radiation  pattern  away  from  the  interface  shows  similarities 
with  an  effective  electric  dipole  radiation  pattern  centered  right  beneath  the  source  at  the 
contrast.  The  TM  mode  amplitudes  decay  rapidly  with  traveled  distance  and  suggests  the 
importance  of  a  Vertical  ElectroSeismic  Profiling  geometry  to  record  the  converted  electro¬ 
magnetic  signal  at  antennas  close  to  the  target  (contrast)  of  interest. 

1.  Introduction 

When  seismic  waves  propagate  through  a  fluid  saturated  sedimentary  material  a  small 
amount  of  relative  fluid-solid  motion  is  induced  (the  motion  of  the  pore  fluid  to  the  solid 
matrix  defines  relative  flow).  The  driving  force  for  the  relative  flow  is  a  combination  of 
pressure  gradients  set  up  by  the  peaks  and  troughs  of  a  compressional  wave  and  by  grain 
accelerations.  The  relative  flow  caused  by  grain  accelerations  can  therefore  be  due  to  either 
compressional  or  shear  waves. 

A  fluid  electrolyte  in  contact  with  a  solid  surface  chemically  adsorbes  the  anions  from  the 
electrolyte  to  the  grain  surfaces  leaving  behind  a  net  excess  of  ions  distributed  near  the  wall. 
This  region  is  known  as  an  electric  double  layer  (Bockris,  1970).  The  diffuse  distribution  of 
mobile  ions,  with  a  higher  concentration  ions  in  the  region  close  to  the  adsorbed  layer  and 
more  and  more  diffuse  towards  the  neutral  electrolyte  are  free  to  move  when  the  fluid  moves. 
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The  seismic  wave  motion  which  generates  the  relative  flow  also  induces  a  ‘streaming’ 
electrical  current  due  to  flow  of  double-layer  ions.  This  induced  streaming  current  acts  as  a 
current  source  in  Maxwell’s  equations.  The  compressional  waves  in  a  homogeneous  porous 
media  contain  an  nonradiating  electric  field  that  is  not  detectable  outside  the  wave  pulse. 
The  peaks  and  troughs  cause  a  charge  separation  setting  up  electric  fields  of  the  Coulombic 
attraction  type,  that  drives  a  conduction  current  that  exactly  balances  the  streaming  current. 
Propagating  electromagnetic  waves  are  not  generated  since  there  is  no  current  imbalance 
within  the  seismic  pulse. 

Divergence-free  S  waves  cannot  cause  charge  to  separate  since  they  do  not  induce  changes 
in  fluid  pressure  but  there  is  relative  fluid-solid  motion  through  equi-voluminal  grain  accel¬ 
erations.  The  induced  streaming  current  sheets  induce  magnetic  fields  that  generate  very 
small,  compared  to  the  electric  field  associated  with  the  P-wave  pulse,  electric  fields  of  the 
induction  type.  When  rotational  waves  propagate  through  homogeneous  porous  medium, 
it  contains  a  measurable  magnetic  field  within  the  pulse  that  is  not  detectable  outside  the 
seismic  pulse.  In  homogeneous  media  a  dipole  antenna  acts  as  a  geophone  measuring  the 
nonradiating  electric  fields  inside  the  compressional  pulse  when  the  pulse  passes  by.  A 
magnetometer  would  act  as  a  shear  wave  selective  geophone  measuring  the  non-radiating 
magnetic  fields  inside  the  rotational  pulse  when  the  pulse  passes  by. 

Seismic  waves  that  traverse  a  contrast  in  electrical  and/or  mechanical  properties  how¬ 
ever  form  a  complex  dynamic  current  imbalance  across  the  interface  generating  radiating 
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electromagnetic  waves  measurable  at  the  earth’s  surface.  Field  data  has  been  recorded 
that  demonstrates  that  when  seismic  sources  which  generate  seismic  waves  that  propagate 
through  the  near  surface  sediment  layers,  convert  at  depth  into  electromagnetic  disturbances 
that  can  be  recorded  (Thompson  and  Gist,  1991,  1993;  Butler  et  al.,  1994;  Mikhailov  and 
Haartsen,  1996). 

In  this  report,  the  macroscopic  equations  controlling  such  behavior  will  be  numeri¬ 
cally  solved  for  the  case  of  seismic  point  sources  in  a  layered  medium.  Only  an  explo¬ 
sive  pointsource  will  be  numerically  modeled,  it  will  be  shown  however  how  to  allow  for 
other  sources  as  well.  A  global  matrix  method  (Chin  et  al.,  1984;  Mai,  1988)  is  employed 
that  solves  simultaneously  for  all  the  macroscopic  electromagnetic  and  poro-elastic  wavefield 
properties.  The  EM  field  converted  from  seismic  waves  traversing  boundaries  are  associ¬ 
ated  with  a  change  in  the  following  properties:  changes  in  elastic  properties  (e.g.,  porosity, 
bulk  and  frame  moduli,  changes  in  fluid  flow  permeability,  changes  in  fluid-chemistry  that 
affect  the  amount  of  double  layer  ions  free  to  move  in  the  double  layer  (e.g.,  bulk  free-ion 
concentration,  pH)). 

A  thick,  relative  to  the  mechanical  wavelengths,  permeable  sand  layer  example  is  pre¬ 
sented  to  show  the  effect  of  a  mechanical  contrast  on  the  mechanical  and  electromagnetic 
wavefield  components.  A  fresh  water/brine  contrast  model  is  used  to  show  the  effect  of 
an  electrical  contrast  on  the  mechanical  and  electromagnetic  wavefield  components.  Both 
models  are  also  used  in  the  calculation  of  converted  magnetic  and  electric  field  amplitudes 
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versus  seismic  source- antenna  offset  at  different  distances  from  the  interface. 

Snapshots  in  time  are  calculated  to  follow  the  conversion  evolution  of  mechanical  waves 
into  electromagnetic  waves.  Xhe  rapid  decay  in  amplitude  with  traveled  distance  of  the 
converted  electromagnetic  signal  driven  by  the  seismic  source  with  a  seismic  center  frequency 
suggests  the  Vertical  ElectroSeismic  Profiling  geometry,  where  antennas  are  positioned  close 
to  the  target  of  interest.  The  last  two  numerical  examples  are  calculations  in  this  VESP 
geometry. 

2.  Macroscopic  Coupled  Electromagnetic  and  Poro-Elastic  Field 
Equations 

Assuming  an  dependence,  Pride  (1994)  derived  the  following  macroscopic  fully  coupled 
mechanical,  electromagnetic  equations  and  constitutive  relations  in  volume  averaged  form 
describing  the  coupled  field  behavior  in  two  phase  porous  medium. 


V-Zb  =  [pbUs  +  Pfm]  +  F 

(1) 

Tg  =  [Kg^  '  Us  +  •  w]  X  -|- 

Vms  +  VmI  -  ^  V  •  UsL 

(2) 

—P  =  CV  •  Ms  +  FIV  ■  w 

(3) 

—iuw  —  L{lo)E_+  ^  [  VP-t-o;  p/Ms  +  /] 

(4) 

J  =  tT(a;)P -1- L(a;)  [— VP  +  a;^P/Ms  + /] 

(5) 
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(6) 


V  X  ^  =  —iujD  +  t/  +  O 


(7) 


D  = 


Co 


—  (k/  -  Ks)  +  Ks 
ftoo 


je; 


(8) 


B  =  noR 


(9) 


were  E_  is  the  electric  field  strength,  H  is  the  magnetic  field  strength,  B  is  the  magnetic 
flux  density,  D  is  the  electric  displacement,  cq  is  the  permittivity  in  free  space,  no  is  the 
permeability  in  free  space,  k/  is  the  relative  permittivity  in  the  fluid,  Kg  is  the  relative  per¬ 
mittivity  in  the  solid,  is  the  bulk  stress,  P  is  the  pressure  in  the  pore  fluid,  Ug  is  the 
displacement  in  the  solid,  w  is  the  relative  fluid-solid  motion  and  I  is  the  identity  matrix.  C 
is  an  external  current  density  source  and,  F  and  /  are  average  bodyforces  densities  acting  on 
the  bulk  and  fluid  phase  respectively.  Equations  (1-4)  characterize  the  mechanical  wavefield 
behavior,  with  Eq  (2)  and  (3)  the  deformation  equations  in  the  fluid  and  solid  phase,  respec¬ 
tively.  Equations  (5-9)  characterize  the  electromagnetic  wavefield  behavior,  with  Eq  (8)  and 
(9)  the  electromagnetic  constitutive  relations.  Equations  (4)  and  (5)  are  the  transport  equa¬ 
tions  through  which  the  mechanical  wavefields  are  coupled  to  the  electromagnetic  wavefields. 
The  coefficients  in  the  flux  equations  relating  the  equations  of  mass,  here  fluid  flow  and  cur¬ 
rent  flow  when  there  are  potential  and  pressure  gradients  are  the  complex  and  frequency 
dependent  conductivity  of  the  two  phase  medium,  <7(0;),  k{u))  is  the  dynamic  permeability, 
L{u))  is  the  complex  electrokinetic  coupling  coeflicient,  and  t]  is  the  fluid  viscosity.  Isotropic 
complex  and  frequency  dependent  expressions  for  coefficients  k{u)  and  L{ui)  are  in 
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Pride  (1994).  The  first  term  in  Eq  (4)  is  the  conduction  current  contribution  and  the  second 
term  is  the  streaming  current  contribution  to  the  total  current  density,  J_.  The  fact  that  the 
cross  terms  have  the  same  coupling  coefficient  L(u))  is  a  statement  of  Onsager  reciprocity. 
The  coefficients  in  the  deformation  equations  are, 


11 

4 

H--G  = 

Kfr  “t“  (j^Kfr  (1  (l))Ks^ 

(10) 

1  + A 

c  = 

K,  +  K,A 

(11) 

1  + A 

M  = 

1  Kf 

1  t  i  A 

(12) 

where  the  parameter  A  is  defined  as, 

A  =  1(1  -  m.  -  KlA  (W) 

(pK; 

The  moduli  Kfr  and  G  are  the  bulk  and  shear  moduli  of  the  framework  of  the  grains,  when 
the  fluid  is  absent.  The  frame  moduli  may  either  be  considered  experimentally  determined 
or  may  be  obtained  from  approximate  theoretical  models  for  specific  pore  grain  geometries. 
C  and  M  are  the  incompressibilities  used  by  Biot  (1962)  and  Pride  (1992),  they  are  complex 
and  frequency  dependent,  allowing  for  losses  in  addition  to  those  associated  with  relative 

flow. 
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3.  Elect roseismic  Wavefield  Description  in  an  Isotropic  Porous  Lay¬ 


ered  Material 


There  are  assumed  to  be  ND  layers.  The  fundamental  input  medium  properties  used  in 
the  numerical  modeling  that  characterize  each  layer  and  determine  all  complex  wavefield 
velocities,  the  complex  and  frequency  dependent  electrical  conductivity  and  electrokinetic 
coupling  coefficients,  are:  porosity,  dc  permeability,  bulk  modulus  of  the  solid  and  fluid 
phase,  the  frame’s  bulk  and  shear  modulus,  the  densities  of  the  solid  and  fluid  phase,  the 
viscosity  of  the  fluid  phase,  the  temperature  of  the  bulk  material,  the  tortuosity,  the  salinity 
of  the  fluid  and  the  relative  permittivities  of  the  solid  and  fluid  phases. 

First,  the  field  equations  describing  the  wave  behavior  in  each  layer  are  determined.  A 
cylindrical  coordinate  system  (r,  ^  ,  z)  with  z  being  depth  has  been  used  since  the  geometry 
has  azimuthal  symmetry.  The  following  definitions  involving  the  horizontal  components  and 
horizontal  derivatives  are  used  (Hudson,  1969;  Kennett,  1983), 


Uy 


d  .  .  d 
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d ,  ,  d 

Tr 
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(15) 
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(19) 


<  E,H  >H= 


1 


^2  Id  f  d' 

^^  =  rd-rVYr 


+ 


dijp' 


(20) 


to  rewrite  the  governing  macroscopic  coupled  electromagnetic  and  poro-elastic  field  equa¬ 


tions.  To  obtain  the  transform-domain  field  equations,  the  governing  equations  controlling 


the  coupled  electromagnetic-seismic,  or  electroseismic  wavefield  propagation  are  first  Fourier 
transformed  to  the  temporal  frequency  domain.  A  Finite  Fourier  Hankel  transform  is  per¬ 
formed  over  the  horizontal  coordinates  r  and  <j)  thus  providing  a  cylindrical  wave  decompo¬ 


sition. 

This  transform  pair  is  defined  as, 

HF  (j))]  =  i>{k,  ^  ^  rJn{kr)  del)  (21) 

HF~^\'^{k,n)]='ip{r,(l>)=  I  dk  k  ^  J„(A:r)e*"^i^(fc,n)  (22) 

'•  ■'  •'o  n=-N 

Note  the  property  =  -k^,  (Watson,  1944).  The  magnitude  N  in  the  summation 

is  determined  by  azimuthal  symmetry  of  the  point  source.  For  spherically  symmetric  sources 
there  is  no  azimuthal  dependence,  then  N  =  0.  If  the  source  can  be  described  using  an 
arbitrarily  directed  force  vector,  then  N  <1  ,  while  if  it  can  be  described  using  a  second 
order  moment  tensor,  then  N  <2  (Kennett,  1983). 

The  final  transform-domain  field  equations  are  obtained  by  eliminating  those  field  compo¬ 
nents  that  are  discontinuous  across  an  interface  between  two  media  with  different  poro-elastic 
and/or  electromagnetic  medium  properties.  The  resulting  differential  equations  are  written 
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in  the  form  of  a  first-order  ordinary  matrix  diflferential  equation.  In  isotropic  media  the  first- 
order  ordinary  differential  equations  decompose  into  two  sets  that  describe  PSVTM  and 
SHTE  coupled  wavefield  propagation.  The  compressional  and  vertical  polarized  rotational 
mechanical  waves  generate  electrical  currents  in  the  PSV  partical  motion  plane  in  homo¬ 
geneous  isotropic  media.  This  couples  to  the  electromagnetic  wavefield  components  with 
TM  polarization.  We  call  this  from  now  on  the  PSVTM  coupled  electroseismic  wavefield 
case.  The  horizontal  polarized  rotational  mechanial  waves  generate  electrical  currents  in 
the  SH  particle  motion  plane  in  isotropic  homogeneous  porous  media.  This  couples  to  the 
electromagnetic  wavefield  components  with  T E  polarization. 

An  incoming  wave  at  an  interface  separating  two  isotropic  porous  materials  will  generate 
four  reflected  and  transmitted  wavetypes  {Pfast,  Psiow,  EV,  TM)  in  each  porous  layer  in  the 
PSVTM  case.  Thus,  eight  boundary  conditions  need  to  be  complied  with,  requiring  the 
phase  matching  of  eight  wavefield  components  across  the  interface.  In  the  SHTE  case,  an 
incoming  wave  at  an  interface  generates  two  reflected  and  transmitted  wavetypes  (SH,  TE). 
Thus,  four  boundary  conditions  need  to  be  complied  with,  requiring  the  phase  matching  of 
four  wavefield  components  across  the  interface.  In  Pride  and  Haartsen  (1995)  the  electroseis¬ 
mic  boundary  conditions  are  derived  from  kinematic  constraints  and  the  conservations  laws 
(linear  momemtum,  conservation  of  mass  and  conservation  of  energy).  The  obtained  bound¬ 
ary  conditions  were  shown  not  to  violate  uniqueness  of  solution  of  an  electroseismic  boundary 
value  problem.  These  derived  boundary  conditions  require  the  following  transform-domain 
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continuity  displacement-stress-EM  wavefield  components  in  layer  m  to  be  continuous  across 
its  interfaces  with  other  layers, 


^  ^  -  iT 

ZJ  ^Z7  ^1)  '^ZZt  -^21  -^1 

(23) 

r,(m;SHTE) 

- 

^U2,f2,Hi,E2] 

(24) 

were  the  horizontal  components  are  defined 

as  follows, 

fiv  ..  '^Vz 

~  ik  ik 

II 

S.; 

II 

(25) 

fi//  -  ‘^Hz 

~  ik  ik 

fi  -Ik  F 

^  ik  ^  ik 

(26) 

The  advantage  of  using  these  horizontal  field  components  becomes  clear  if  one  wants  to 
solve  a  line  source  parallel  with  the  stratification  instead  of  a  pointsource.  Working  in 
the  {x,  y,  z)  cartesian  coordinates,  assuming  the  line  source  to  be  parallel  with  the  y  axis 
and  taking  a  Fourier  transform  of  the  equations  of  motion  with  respect  to  x,  a  plane  wave 
decomposition  is  now  employed,  results  in  exactly  the  same  equations  as  are  obtained  when 
the  following  identification  of  the  horizontal  components  in  the  displacement-stress-EA/ 
vectors  are  made:  ui  =  fix,  fia  =  fiyi  ==  '^2  =  ^yz^  Hu  =  Hy,  Hy  —  Hx-,  Eh  —  Ey  and 

Ey  -  Ex-  Thus,  using  this  mapping,  it  follows  that  the  solution  of  the  first-order  differential 
equations  applies  directly  to  the  line  source  problem  as  well.  Additionally,  the  cylindrical 
wave  reflection/transmission  coefficients  to  be  obtained  are  necessarily  identical  to  the  plane 
wave  coefficients. 


110 


The  PSVTM  continuity  vector,  Eq  (25)  is  related  to  the  medium  at  all  depths  z  through 
an  eight  by  eight  system  matrix  The  SHTE  continuity  vector,  Eq  (26)  is  related 

to  the  medium  at  all  depths  2:  through  a  four  by  four  system  matrix  The  system 

matrix  components  depend  on  frequency  and  complex  frequency  dependent  poro-elastic, 
electric  and  magnetic  medium  properties,  the  horizontal  slowness  p  (note  k  =  uip  )  and  the 
transport  and  electrokinetic  coupling  coefficients.  After  Hankel  transforming  Eq  (6-3)  and 
some  algebra  the  following  first-order  linear  differential  matrix  equations  which  describe  the 
PSVTM  case  are  obtained. 


d  n(PSVTM) 


PSVTM) 


A{osmosis;PSVTM) 


A  (electr akinetic; PSVT M) 


A(Maxwell;PSVT  M) 


r.{PSVTM) 


(27) 


The  first  order  differential  matrix  equation  describing  the  SHTE  case  in  the  transform 
domain  is. 
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where  the  6  by  6  PSVTM  submatrix  and 
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relate  the  mechanical  stress-displacement  wavefields  with  depth. 

Where  the  2  by  2  PSVTM  submatrix  and  the  2  by  2  SHTE  submatrix 


AMaxwell-,PSVTM)  _ 

^IJ  ~  §33’ 


AMaxwell\SHTE)  _ 

Aj 


0  -io,  (l  -  g) 


—iuH  0 

relate  the  electric  and  magnetic  wavefield  components  with  depth. 

Where  the  6  by  2  PSVTM  submatrix  and  the  2  by  2  SHTE  submatrix 


j  (electrokin€tic\PSVT M) 


t(electrokinetic‘,SHTE)  _ 


u^Lpf  0 
0  0 


(30) 


(31) 


relate  the  electric  and  magnetic  wavefield  components,  converted  from  the  mechanical  wave- 


fields,  with  depth. 

And  where  the  2  by  6  PSVTM  submatrix  and  the  2  by  2  SHTE  submatrix 


1  (osmosis;P  SVT  M) 


il3 


=23 


j^{osmosis;btH  _  q  (32) 

relate  the  mechanical  wavefield  components,  converted  from  electric  and  magnetic  wavefields, 
with  depth. 

The  following  matrices  have  been  used, 
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(34) 


B  -  2p^G{l  +  a) 

0 

0 

0 

0 

§21  = 

0 

Pb 

Pf 

’  §23 

0 

0 

0 

Pf 

f 

0 

P^€-peL^ 

Lpf  0  0 

0  0  -2^ 
€-PeL^ 


’  §32 


ILOe 


0  0  L 
0  0  0 


The  following  coefficients  have  been  used  in  the  matrices  above. 
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4.  The  Global  Matrix  Method 

To  obtain  the  wavefields  in  each  layer  m,  the  following  continuity  conditions  at  each  boundary 
are  used, 

lim  —  lim  =0  s  (39) 

ziZm  zJZm 

At  the  source  level  z  =  Zs,  where  the  source  is  positioned,  Bj  jumps  by  a  finite  amount.  At 
the  source  level, 

lim  -  lim  B^  =  Sj  (40) 

zl,Zg  Z^Zg 
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with  Sj  the  source  vector.  First  a  linear  transformation  on  each  of  the  ND  field  vectors 
is  carried  out.  Through  it,  a  field-vector-formalism  is  obtained  in  which  a  decomposition  of 
into  up  and  downgoing  fields  is  manifest.  Let  eight  by  one  vector 

and  a  four  by  one  vector  that  is  related  to  by  the  linear  transformation, 

(41) 

where  matrix  is  subject  to  a  convenient  choice.  From  the  theory  of  matrices  Dj^j^ 
is  taken  the  local  eigenvector  matrix  for  the  system  matrices  defined  in  Eq  (27)  and 
(28).  In  the  PSVTM  case  the  eigencolumn  matrix,  for  each  layer  m,  of  the  system  matrix 
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^P,Q  ^^^P,Q 


(42) 

(43) 


j,r{m, PSVTM-,-)  ,r{m,PSVTM-,+) 

[  ^^P,Q  ^^P,Q  J 

where  d  contains  the  eight  by  one  eigenvectors  with  r  G  Pfast,Psiow,SV,TM,  the  up 

and  downgoing  wavefields,  shown  in  Eq  (A-4)  and  (A-7)  in  Appendix  A.  Mpq  ^  and 

j\j(^^P^VTM;±)  jnatrices. 

In  the  SHTE  case  the  eigencolumn  matrix,  for  each  layer  m  of  the  system  matrix 
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where  d  contains  the  four  by  one  eigenvectors  with  r  G  SH,TE,  the  up  and  downgoing 


wavefields,  shown  in  Eq  (A-6)  in  Appendix  A. 

Substituting  Eq  (41)  into  the  first-order  differential  matrix  equation,  Eq  (27)  and  (28), 


and  premultiplying  the  equation  with  ^  yields, 


—  i  Aj  j  Uj  j^  -  Um,I  ■ 


In  a  uniform  medium  the  eigenvector  components  in  are  independent  of  z,  with  the 


result  that  ^ vanishes.  Since  Dj^n  is  taken  the  local  eigenvector  matrix  for  the 


system  matrices  we  obtain, 
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where  A^*jv  is  a  diagonal  matrix  of  the  eigenvalues  of  the  system  matrices 
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The  vector  is  governed  by  the  following  differential  equation, 


with  solution. 
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We  define  the  following  wave-propagator  matrix, 


expliuq^’^'  '^\z  -  Zm)] 


with  the  matrices  defined  as, 


0 

0  4”^-‘ 


(51) 


=  diag  (exp[iuq^’^'’^^^Az],exp[iuq^"^'’^‘'^Az], 
exp[iu}q^'^'^^^  Az],  exp[iuq^"^’’^^^  Azfj 

=  diag  (exp[iu)q^”'''^^'^  Az],  exp[iuq^'^'’'^^'> Az]^ 


(52) 

(53) 


with 

vertical  slownesses  for  the  compressional  fast  wave  mode,  the  compressional  slow  wave  mode. 


the  vertical  polarized  shear  wave  mode,  the  horizontally  polarized  shear  wave  mode,  the 


transverse  magnetic  EM  mode  and  the  transverse  electric  EM  mode.  The  choice  of  bi  anch 
cuts  is  defined  according  to  the  computational  methods  used  to  perform  the  inverse  Fourier 
and  Hankel  transforms;  >  0,  >  0  and  >  0. 

The  field  vector  wP{z)  in  layer  m  is  related  to  field  vector  1T^’”'*'^^(2:)  in  layer  m  +  1 
through  the  continuity  conditions  shown  in  Eq  (39)  and  (40).  The  interface  continuity 
conditions  for  m  ^  s  after  applying  the  linear  transformation  are, 
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^P,Q 
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At  z  =  zyv  the  interface  conditions  are, 


^P,Q  ^P.Q  t  ^P,Q  P 


^(NA-yND;-)^0 


To  include  a  free  surface,  the  free  surface  conditions  for  the  mechanical  wavefields  and  the 
ordinary  electromagnetic  boundary  conditions  at  the  porous  medium-air  interface  need  to 
be  complied  with.  The  mechanical  normal  and  tangential  stresses  and  pressure  in  the  two 
phase  medium  have  to  vanish  at  the  free  surface,  while  the  tangential  electric  and  magnetic 
fields  obey  the  porous  medium-air  boundary  conditions.  When  we  neglect  the  osmosis  effect 
at  the  free  surface,  the  up  and  downgoing  fast,  slow  and,  shear  wave  amplitudes  can  be 
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related  through  the  vanishing  of  the  stress  and  fluid  traction  boundary  conditions.  The  up 
and  downgoing  electric  and  magnetic  wavefields  at  the  free  surface  are  related  through  the 
ordinary  electromagnetic  boundary  conditions.  All  downgoing  wave  type  amplitudes  can  be 
related  to  the  upgoing  wave  amplitudes  and  satisfy  all  free  surface  conditions.  The  reflected 
wavefield  from  the  free  surface  are  phase  delayed  to  the  interface  below  and  included  in  the 
phase  matching  across  this  boundary. 

In  this  analysis  the  following  matrices  are  used.  The  3x3  matrix,  which  contains 
the  mechanical  eigenvector  components  that  relate  the  mechanical  wavefield  amplitudes  to 
the  normal  and  shear  bulk  stresses  and  pressure  in  the  pore  fluid.  The  2x3  matrix, 
which  contains  the  mechanical  eigenvector  components  that  relate  the  mechanical  wavefield 
amplitudes  to  the  tangential  magnetical  and  electrical  wavefields.  The  2x2  matrix,  1^2^ 
which  contains  the  electromagnetic  eigenvector  components  that  relate  the  EM  wavefield 
amplitudes  to  the  tangential  magnetical  and  electrical  wavefields  in  the  porous  medium. 
The  2x1  matrix,  W2i'~^  which  contains  the  electromagnetic  eigenvector  components  that 
relate  the  EM  wavefield  amplitudes  to  the  tangential  magnetical  and  electrical  wavefields  in 
the  upper  halfspace  (air).  Where  the  superscript  1  denotes  the  first  layer,  0  the  halfspace 
above  the  free  surface,  and  T  denotes  eigenvector  components  relating  to  the  down  or  up 
going  wavefields.  Using  the  eigencolumn  components  listed  in  Appendix  A,  the  following 
matrices  are  constructed,  necessary  to  express  the  reflected  wavefield  amplitudes  from  the 
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free  surface  in  terms  of  upgoing  wave  amplitudes. 
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(56) 


The  up  and  downgoing  mechanical  wavefield  amplitudes  in  the  first  layer,  are  the  fast  wave 
amplitudes,  the  shear  wave  amplitudes,  and  the  slow  wave  amplitudes, 

which  can  be  related  through  the  vanishing  stress  and  fluid  traction  boundary  conditions  at 
the  free  surface. 


jn-) 

7^1-3 


c(i;+)Ji;+) 
*^33  Fi-3 


(67) 
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The  electromagnetic  boundary  conditions  are  the  usual  continuity  of  the  magnetic,  and 
electric,  Er,  components  of  the  electromagnetic  T M  mode.  The  radiation  in  the  space  above 
the  free  surface  requires  the  condition  of  no  incoming  EM  waves  from  ^  —  oo,  implying 

=  0,  downgoing  TM  wave  amplitude  is  zero.  The  EM  boundary  conditions  at  the  free 
surface  are  expressed  as  follows. 


^'2.1  P\-3  +1^23  Pl-3  +*'2 


22 


=  w. 


(0;+) 


21 


(58) 


Using  Eq  (57),  the  following  relation  between  the  up  and  downgoing  EM  wave  amplitudes 
propagating  in  the  first  layer  can  be  obtained. 


=  [n«] (59) 

J 

-  + 41" VS>] ) 

Eliminating  the  amplitude  ,  the  downgoing  TM  amplitude  in  the  first  layer  can  be 

expressed  in  terms  of  the  upgoing  T M  and  mechanical  wavefield  amplitudes. 


Pi  '  ^ 


- ^y-  [e  +  H(2,  +  E(2,  2)pS‘-+^  +  E(2, 3)p^'=+^) 

det{y22  ) 
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With, 
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Equations  (57)  and  (60)  relate  the  upgoing  wave  amplitudes  to  the  downgoing  wave  ampli¬ 
tudes  satisfying  the  free  surface  conditions.  The  downgoing  wave  fields  reflected  from  the 
free  surface  are  phase  delayed  to  the  next  interface  and  included  in  its  boundary  condition. 

Collecting  all  interface  condition  equations  in  one  matrix,  a  (ND-l)  x  8  by  (ND-1)  x  8 
dimensional  global  matrix  is  obtained  in  the  PSVTM  case  and  a  {ND—l)x4  by  (7VD— l)x4 
dimensional  global  matrix  in  the  SHTE  case.  To  obtain  the  up  and  downgoing  amplitudes 
in  all  layers  simultaneously  in  both  PSVTM  and  SHTE  cases,  the  global  matrix  needs 
to  be  inverted.  The  obtained  global  matrix  is  identical  in  shape  to  the  matrix  obtained 
for  a  radiating  upper  halfspace.  Inverting  the  global  matrix  yields  only  the  upgoing  wave 
amplitudes  in  the  first  layer.  The  downgoing  wave  fields  are  calculated  using  Eq  (57)  and 
(60).  The  physical  wavefields  at  the  receivers  are  obtained  by  substituting  the  inverted  up 
and  downgoing  amplitudes  in  each  layer  into  Eq  (41)  and  propagate  them  to  the  receiver 
location  in  the  layer. 

The  global  matrix  formulation  was  developed  by  Chin  et  al.  (1984).  The  global  matrix 
formulation  has  several  important  features:  (1)  Multiple  sources  can  be  treated,  the  produced 
wavefields  are  simply  superposed;  (2)  The  up  and  down  going  wavefield  coefficients  are 
simultaneously  solved  for  in  all  layers,  therefore  the  wavefield  is  known  to  all  receivers  at  all 
places  in  the  medium;  (3)  Time  stability  problems  do  not  occur  because  decaying  exponents 
are  needed  only.  The  Thomson-Haskell  propagator  matrix  approach  (Schmidt  and  Tango, 
1986)  does  not  have  these  advantageous  features. 
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The  global  matrix  system  that  needs  to  be  solved  is, 
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Chin  et  al.  (1984)  show  that  the  Thomson-Haskell  method  (Thomson,  1950;  Haskell,  1953) 
can  be  viewed  as  treating  the  coefficient  matrix  of  the  linear  system  in  a  block  diagonal  form 
and  applying  a  forward  marching  algorithm  for  its  solution.  Chin  s  formulation  is  a  shooting 
method  for  a  two  point  boundary  value  problem.  Gaussian  elimination  with  partial  pivoting 
yields  stable  solutions.  To  avoid  time  stability  problems  (exploding  exponentials)  the  Wm 
vector  is  multiplied  with  a  diagonal  phase  matrix  Fu  that  cancels  the  potentially  unstable 
upgoing  wave  phase  factors. 


Fu  =  diag  Ip,q,  Ep)^, ...,  E^j^q\Ip,q\ 


(63) 
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where  /p,q  is  the  identity  matrix.  If  we  define  Wm  =  Fmn^n,  can  be  shown  to  satisfy, 
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All  terms  that  contained  Xz] 

are  eliminated.  These  exponentials  cause  numerical  instability  when  the  wavefield  goes 
evanescent  (Dunkin,  1965;  Kennett,  1983;  Schmidt  and  Tango,  1986). 

5.  Point  Forces  in  a  Poro- Elastic  Stratified  Medium 

To  model  a  general  source,  the  moment  tensor  representing  a  body  force  equivalent  in  a 
poro-elastic  medium  has  been  used.  Since  the  effect  of  the  source  is  an  internal  process 
within  a  volume  V,  its  total  momentum  and  total  angular  momentum  must  be  conserved 
(Aki  and  Richards,  1980).  Therefore,  the  total  force  and  momentum  about  any  fixed  point 
must  be  zero  (i.e.,  fy  h{r)dV{r)  =  0  and  fy  h(r)  x  (r-r^)dV(r)  =  0,  with  h  a  body  force  and 
Zg  a  fixed  point  in  space).  The  above  conditions  imply  the  following  body  force  equivalent, 
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hj  =  -Ls)  Mjk  =  Mkj. 

The  following  equation  has  to  be  evaluated  when  a  cylindrical  coordinate  system  is  used, 


5.1.  Mechanical  Point  Sources 


First  an  explosion  source,  vertical  line  source  and  vertical  dipole  source  uniformly  acting  on 
the  solid  frame  and  pore  fluid  are  derived.  These  sources  have  due  to  symmetry  no  angular 
dependence.  Then  a  couple  source  acting  in  the  2  direction  with  a  spacing  between  the  forces 
in  the  r  direction,  having  angular  dependence  is  outlined.  Sources  uniformly  acting  on  the 
two  phase  medium  are  given  by, 

h  =  mV6{r  — Zs)  =  mV  -^6{(j))6{z  —  Zs)  (65) 

with  H  a  source  acting  on  the  bulk  material  and  h  a  source  acting  on  the  pore  fluid  phase 
of  the  porous  medium. 

In  component  form  the  explosion  source  is  given  by. 


Vi 
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1 - 1 

II 

(66) 

(67) 

K  =  m  ^  ^S{(j))  [(5(^-2,)] 

r  oz 

(68) 

124 


The  horizontal  component  of  the  explosion  source  using  definition  14  reads, 


hv  — 


d  .  .  d 


=  m6(z  — 


6{r)6{(l>) 


(69) 


The  horizontal  and  vertical  component  of  the  explosion  source  after  Hankel  transformation 
and  mapping  on  a  plane  wave  description  read, 

hi  =  =  ikm6{z  —  Zs)  (70) 

tup 

hz  =  m^b{z  -  Zs)  (71) 

az 

The  vertical  line  source  components  are  given  by  70  and  hz  =  0  and  the  vertical  dipole  force 
is  given  by  hi  =0  and  71. 

A  vertical  couple  source  in  the  z  direction  with  a  spacing  in  between  the  forces  in  the  r 
direction  is  given  by, 


d 

hr  =  h^  =  0,  hz  =  m— 


6{r) 


6^— 6(2  -  2s)- 


(72) 


The  vertical  component  of  the  vertical  couple  source  after  Hankel  transformation  (Sneddon, 
1951),  and  mapping  on  a  plane  wave  description  read. 


hi  =0,  hz  =  ±—8{z  —  Zs)  (73) 

oz 

The  effect  of  a  general  mechanical  point  source  in  a  stratified  porous  medium  is  accommo¬ 
dated  by  specifying  a  jump  in  the  displacement  stress  vector  Bj,  across  a  horizontal  plane 
containing  the  source  (Hudson,  1969;  Kennett,  1983). 
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To  obtain  this  jump  vector,  the  governing  poro-elastic  constitutive  Eq  (2)  and  (3),  and 
the  Fourier  transformed  equations  of  motion  with  a  point  source  uniformly  acting  on  the 
solid  frame  and  pore  fluid  in  a  poro-elastic  medium, 

-u)^[pBn  +  Pfm]=^-Z  +  H  (74) 

-uj^  [Pf2L  +  Pew]  =  -  VP  +  h  (75) 

are  Hankel  transformed  and  manipulated  into  a  set  of  equations  that  have  all  derivatives 
with  respect  to  z  on  the  left  hand  side,  and  a  right  hand  side  identical  to  the  mechanical 
submatrix  of  the  electroseismic  system  matrix  29.  The  first  set  describes  the  change  in  2:  of 
the  field  quantities  [ui,  Uz,w„  fj,  S] ,  where  the  mapping  given  in  Eq  (25),  (26)  is  applied 
to  describe  the  Pj  —  Pg  —  SV  case, 


where  is  defined  in  Eq  (29).  To  include  some  form  of  source  in  the  stratification 

we  must  solve  the  following  inhomogeneous  equation, 

B{z)  =  F{z)  (77) 
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First  we  multiply  Eq  (77)  through  with  a  propagator  matrix  R  ^{z,Zs)  from  the  left  hand 
side  to  obtain  an  integration  factor. 


£  ^(2,  2s) 


dz 


B{^) 


B(z)=g(z.,z)F(z) 


(78) 


where  we  have  used  the  relationship  £(2,,  2)  =  £  ^(2, 2s). 

Then  we  integrate  with  respect  to  2  and  multiply  with  R{z,Zs),  we  obtain, 


B{z)  =  £(2,  Zs)B{zs)  +  r  dC  £(2,  C)£(C) 


(79) 


A  point  source  can  be  manipulated  into  some  dipolar  contribution  (Kennett,  1983). 
This  description  is  given  by, 

F,(z,)  =  Fi'h{z  -  z.)  +  -  z.) 


(80) 


with  F/**  and  F,^*  the  dipolar  contributions.  Substituting  Eq  (80)  into  Eq  (79),  we  obtain 

/'dCfifeO 

Jzs 

H{z  -  Zs)Eiz,  2s) 


fI"6{z  -  z.)  +  fP^S(z  -  z.) 
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E{z.,z.)F<'^  -  ^Eiz.X)  £'"’ 


(81) 


with  H{z)  the  Heaviside  step  function. 

Using  the  relationships,  £(2s,C)  =  Er^{z,Zs)  and  ^s)  =  —E  ^{z,  Zs)A{p,  z),  we 

obtain  a  displacement-stress  jump  vector  condition  across  a  source  plane. 

Si{zs)  =  $(uj)  ^Bi{zf)  -  5/(27)]  =  [£/*^  +  2s)£/  (82) 


with  0(a;),  the  Fourier  transform  of  the  source  signature  (t){t). 
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The  final  displacement-stress  jump  representations  of  some  mechanical  point  sources  in 
a  poro-elastic  isotropic  medium  are, 

Explosive  point  source; 


(83) 
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Vertical  couple  source: 


0,  ±ijjpnS^  — ,  (h-C-  ,  ±iu^p'^m  (/?  -  a) ,  0, 0 

A  A  V  Pe  ) 


(85) 


Vertical  dipole  source: 


Mui) 


(86) 


5.2.  Electromagnetic  Point  Sources 


The  effect  of  an  electrical  point  source  in  a  stratified  porous  medium  is  accommodated  by 
specifying  a  jump  in  the  electromagnetic  field  component  vector  ^2)  ^i]  across  a  horizontal 
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plane  containing  the  source.  To  obtain  this  jump  vector,  Maxwell’s  equations  with  an  exter¬ 
nal  current  source  C  in  Fourier  transform  domain  are  Hankel  transformed  and  manipulated 


A  dipole  current  source  produces  an  electromagnetic  field  that  is  expected  to  generate  seismic 
waves  when  it  traverses  an  interface  at  depth.  In  Pride  and  Haartsen  (1995)  it  is  shown, 
using  a  reciprocity  argument,  that  in  order  to  detect  a  seismic  displacement  of  the  order  of 
10~^  m  a  10  Cm  dipole  moment  is  needed  to  produce  an  electric  field  on  the  order  of  10® 
V/m  at  a  distance  10  m  from  the  source.  Electromagnetic  waves  do  generate  mechanical 
wavefields  at  depth  that  are  possibly  measurabe  at  the  earth’s  surface.  This  very  small 
conversion  from  electromagnetic  waves  into  mechanical  waves  can  be  numerically  simulated. 


6.  Transformation  Back  to  the  Space  Time  Domain 

In  this  section  the  back  transformation  to  the  space-time  domain  of  the  displacement-stress- 
EM  vector  is  performed.  An  inverse  Fourier  transform  is  applied  to  go  back  to  the  time 


129 


domain.  An  inverse  Hankel  transform  is  applied  to  obtain  the  3D  spatial  dependence  of  the 
displacements,  stresses,  electric  and  magnetic  fields.  The  horizontal  components  of  the  dis¬ 
placements,  stresses,  electric  and  magnetic  fields  require  additional  integration  over  r  and  0 
to  obtain,  u,.,  r^zi  'T’rzi  Hti  Htf,  and  Er,  due  to  the  definitions  of  uy ,  uh ,  Ty/^,  tjjzi  Hy ,  Hh 


and  Ey,  Eh- 

Uz{u,r,(l),z)=  [  kdk 

*'0  n=-N 

Identical  relations  exist  for  Wz,Tzz  and  S  =  —P.  The  horizontal  components  may  be  recov¬ 
ered  using  the  following  equations  (Kennett,  1983), 


■M, 


POO  Jl 

.{u},r,4>,z)  =  /  kdk  —Jn{kr)u2{to,k,n,z) 

-  iJn{kr)ui{u),k,n,z)]e‘'''^ 


(90) 


POO  Ti 

Ju,r,(j),z)  =  /  kdk  Y 

Jo  n=-N 

-f-  iJ'^{kr)u2{co,k,n,z)]e'’^‘^ 


(91) 


Identical  relations  exist  for  r^z,  t<^z  terms  of  fi  and  f2,  for  E-ry  Eff,  in  terms  of  Ei  and  E2  and 
Hr,  in  terms  of  ^1  and  ^2-  The  above  representations  may  be  regarded  as  a  superposition 
of  cylindrical  waves  whose  order  dictates  the  nature  of  their  azimuthal  modulation.  At 
each  frequency  and  angular  order  the  radial  contribution  is  obtained  by  superposing  all 
horizontal  wavenumbers  k.  This  corresponds  to  including  all  propagating  waves  at  the  level 
within  the  stratification,  from  vertically  to  purely  horizontal  traveling  waves  including  the 
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evanescent  waves.  At  any  particular  distance  r  the  relative  contributions  of  the  wavenumbers 
are  imposed  by  the  radial  phase  functions  Jn{kr). 

The  integrals  (89),  (90),  and  (91)  are  evaluated  by  the  discrete  summation  over  many 
wavenumbers  or  the  so  called  Discrete  Wavenumber  Method  (Bouchon  and  Aki,  1977). 

The  discretization  of  the  radial  wavenumber  k  in  cylindrical  coordinates  introduces  pe¬ 
riodicity  into  the  source  distribution.  The  original  single  source  problem  changes  after  dis¬ 
cretization  into  periodic  concentric  sources  around  the  original  source.  The  periodicity  of 
these  sources  or  the  distance  between  two  adjacent  circular  sources,  L,  is  related  to  the 

discretization  interval  of  the  wavenumber,  Ak,  by  the  sampling  relation, 

27r 


L  = 


Ak 


(92) 


L  and  therefore  Ak  is  determined  by  assuming  a  receiver  located  at  Xr  =  (ro,  zq)  and  a  source 
at  Xg  =  (0,  Zg)  on  the  symmetry  axes  of  the  medium  configuration.  Given  the  time  window 
to  record  radiated  waves  from  0  to  tman  ’pseudo’  waves  radiated  from  the  periodic  sources 
are  not  allowed  to  enter  this  time  window.  This  requirement  is. 


'\J (Z/  “1“  (*^0  ^  '^fastestir, 


(93) 


or 


i  >  '0  +  -  (20  -  2.)^ 


(94) 


The  sampling  equation  becomes  now, 

Ak  < 


27r 


To  +  \lv)aslest^liax  “  (^0  “  Zgf 


(95) 
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The  above  equation  is  the  criterion  for  choosing  the  sample  rate  for  the  discrete  wavenumber 
summation. 

To  perform  the  summation  the  singularities  of  the  integrands  are  moved  further  from  the 
real  k  axes.  This  is  done  by  adding  a  small  imaginary  part  to  the  real  frequency  (analytic 
continuation  into  the  complex  plane),  i.e.,  uj  =  ujr  +  iui  with  uj  >  0.  If  e*"*  was  used 
in  the  time  to  frequency  transformation  there  would  be  a  minus  sign  for  the  imaginary 
part  of  the  frequency.  The  effect  of  complex  frequency  moves  the  singularities  into  the  first 
and  third  quadrant  of  the  complex  A:-plane.  The  use  of  complex  frequency  has  the  effect 
of  smoothing  the  spectrum  and  enhancing  the  first  motions  relative  to  later  arrivals.  This 
effective  attenuation  is  used  to  minimize  the  influence  of  the  neighboring  fictitious  sources 
introduced  by  discrete  k.  The  effect  of  the  imaginary  part  of  the  frequency  can  be  removed 
from  the  final  time  domain  solution  by  inverse  complex  Fourier  transform  with  the  complex 
frequency  with  the  same  imaginary  part  used  in  the  argument  of  the  exponential  function 
in  the  Fourier  transform.  The  magnitude  of  the  imaginary  part  is  usually  chosen  to  be, 

a;;  =  —  (96) 

imax 

Larger  a;/  increases  the  attenuation  for  later  arrivals,  but  also  magnifies  the  numerical  noise 
for  late  times.  If  u>i  is  chosen  too  small,  the  attenuation  may  not  be  large  enough  to  damp 
out  the  arrivals  from  the  fictitious  sources. 
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7.  Numerical  Electroseismogram  Examples 


Four  different  numerical  models  are  calculated  in  the  PSVTM  wavefield  picture.  The  first 
example  is  a  100  in  thick  porous  sand  layer  sandwiched  in  between  two  identical  halfspaces 
that  are  less  porous  than  the  layer.  The  fluid  chemistry  is  the  same  in  all  three  layers. 
The  source/ receiver  medium  configuration  is  depicted  in  Figure  1.  This  example  is  meant  to 
study  electroseismic  conversions  as  a  result  of  a  change  in  mechanical  properties.  Mechanical 
displacement  seismograms  and  electroseismograms  (TM  mode  components)  are  calculated 
and  displayed  in  Figure  2.  The  converted  magnetic  and  electric  field  amplitude  behavior  as 
function  of  antenna  position  are  determined  at  three  different  distances  from  the  mechanical 
contrast  as  shown  in  Figure  3.  In  Figures  4  and  5,  four  time  evolution  snapshots  are  shown 
of  the  mechanical  wavefront  traversing  an  interface  and  the  conversion  of  the  mechanical 
wavefield  into  the  electromagnetic  TM  mode  components. 

The  second  example  studies  the  effect  of  a  change  in  fluid  salinity,  which  only  affects 
mainly  the  medium’s  electrical  properties,  on  the  conversion  to  electromagnetic  waves.  The 
source/receiver  medium  configuration  is  depicted  in  Figure  6.  The  mechanical  displacement 
seismograms  and  electroseismograms  (TM  components)  are  calculated  and  displayed  in  Fig¬ 
ure  7.  The  converted  magnetic  and  electric  field  amplitude  behavior  as  function  of  antenna 
position  at  three  different  distances  from  the  electrical  contrast  are  shown  in  Figure  8.  The 
four  snapshots,  Figures  9  and  10  in  time  are  calculated  to  trace  in  time  the  converted  T M 
wavefield  pattern  generated  by  a  mechanical  wavefront  traversing  the  electrical  contrast. 
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The  last  two  examples  show  the  electroseismic  effect  in  a  Vertical  Electroseismic  Profiling 
(VESP)  setting.  The  seismic  source  is  located  in  the  upper  halfspace  and  the  recording 
geophone/antennas  are  positioned  vertically  crossing  two  mechanical  contrasts  in  the  first 
VESP  example  (see  Figure  11)  and  an  additional  electrical  contrast  in  the  second  VESP 
example  (see  Figure  13). 

The  modeled  mechanical  source  in  all  examples  is  a  mechanical  explosion  source.  To 
model  true  amplitudes,  the  diagonal  elements  in  the  seismic  moment  matrix,  equation  64, 
need  to  be  replaced  by  a  realistic  value.  Using  the  relation  m  =  ^^/CgA©,  with  o  the 
radius  of  the  of  the  nonlinear  deformation  zone  around  the  source,  Ka  a.  coefficient  from 
the  deformation  equation  given  in  equation  10,  A0  the  fractional  change  in  volume  and 
E  =  l2^/i'g.(A©)2  =  e[4.7  X  10^ J/Kg]C,  with  E  the  total  energy,  C  the  source  weight  in 
kilograms  and  e  an  efficiency  factor  expressing  the  fraction  energy  converted  into  sound,  the 
following  expression  for  m  is  obtained. 

m  =  ^^A'g(4.7  X  WJIK9\C€  (97) 

In  the  numerical  calculation  m  =  4.4  x  10^7  is  used,  corresponding  to  C  =  iKg,  a  =  Im, 
A'g  =  lO^Pa  and  e  =  10”^. 

7.1.  The  Electroseismic  Conversion  at  Mechanical  Contrasts 

An  explosion  source  is  positioned  100  m  above  the  first  mechanical  contrast  in  the  upper 
halfspace.  Fifteen  receivers  are  positioned  symmetrically  in  a  straight  horizontal  line  at 
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both  sides  of  the  source  95  m  above  the  interface.  The  receiver  spacing  is  10  m.  The 
medium  parameters  of  the  two  halfspaces  and  the  100  m  thick  sand  layer  are  given  in  Table 
I.  Based  on  these  medium  properties  the  mechanical  fast  wave,  slow  wave,  shear  wave  and 
electromagnetic  TM  wave  velocities  are  calculated.  The  complex  velocities  and  calculated 
bulk  conductivities  are  given  in  Table  II.  A  slow  wave  and  TM  wave  velocity  range  is  given 
since  these  two  wave  phenomena  are  diffusive. 

When  an  explosion  source  is  set  off  only  P  waves  are  generated  in  an  isotropic  poro- 
elastic  medium.  Numerically  the  8  x  1  displacement-stress-EM  wavefield  component  vector  is 
calculated  at  each  geophone/antenna  position.  In  Figure  2  only  the  mechanical  displacement 
seismograms  {u^  and  Ur  components)  and  the  electroseismograms  (T M  mode  components) 
are  shown.  All  plots  are  seismogram  scaled.  The  amplitudes  of  the  first  40  msec  in  the 
Er  and  electroseismograms  are  multiplied  by  a  factor  50  to  enhance  the  electroseismic 
converted  field  generated  at  the  first  interface. 

At  the  mechanical  contrast  the  P  wave  reflects  as  a  F  wave,  a  converted  SV  wave,  and 
a  converted  TM  wave.  The  direct  wave  is  omitted  in  all  seismograms.  The  mechanical 
seismograms  show  the  PP  reflection  and  PSV  conversion  generated  at  the  top  and  bottom 
of  the  thick  sand  layer.  The  converted  TM  wavefield  components  show  up  at  all  antennas 
at  approximately  half  the  two  way  P  wave  travel  time  for  normal  incidence  reflection.  Since 
the  TM  wavefield  velocity  is  at  least  two  orders  of  magnitude  faster  than  the  F  wavefield 
(see  Table  II),  the  traveltime  spent  by  the  TM  waves  in  traveling  upward  to  the  antennas 
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is  negligible  and  the  total  two  way  traveltime  appears  to  coincide  with  the  oneway  trav¬ 
eltime  of  the  incoming  P  wave  at  normal  incidence  reflection.  The  hyperbolas  arriving  at 
a  later  time  in  the  Er  and  electroseismograms  have  the  same  travel  time  as  the  PP 
reflection  (compare  the  mechanical  seismogram  with  the  Er  electroseismogram  in  Figure  2 
and  the  same  traveltime  as  the  PSV  conversion  (compare  mechanical  seismogram  with  the 

electroseismogram),  respectively. 

In  the  first  case  the  compressional  waves  traveling  through  homogeneous  porous  medium 
cause  pressure  gradients  which  cause  charge  to  separate.  This  induces  within  the  seismic 
pulse  a  system  of  electric  fields  that  travel  with  the  compressional  wave  speed.  As  explained 
before,  the  seimic  pulse  doesn’t  radiate  electromagnetic  fields  away  from  the  pulse.  Therefore 
when  the  reflected  P  wave  pulse  passes  the  antennas  an  electric  field  is  registered  inside  the 
P  wave  pulse. 

In  the  second  case  the  vertical  polarized  rotational  mechanical  waves  traveling  through 
homogeneous  porous  medium  cause  grain  accelerations,  setting  up  current  sheets.  This 
induces  within  the  seismic  pulse  magnetic  fields  that  travel  with  rotational  mechanical  wave 
speed.  Therefore  when  the  converted  SV  wave  pulse  passes  the  antennas  a  magnetic  field  is 
recorded  inside  the  SV  wave  pulse.  The  electric  and  magnetic  field  strengths  in  the  seismic 
pulse  are  in  this  case  larger  than  the  converted  TM  wavefield  components  strengths.  The  Er 
hyperbolas  arriving  after  the  electroseismic  conversion  are  associated  with  the  electric  fields 
inside  the  P  wave  pulse  traveling  with  fast  wave  velocity.  While  the  hyperbolas  arriving 
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after  the  electroseismic  conversion  are  associated  with  the  magnetic  fields  inside  the  SV  wave 
pulse  traveling  with  S  wavespeed  velocity.  The  TM  conversions  from  the  bottom  mechanical 
contrast  have  amplitudes  which  are  too  small  to  be  identified  in  the  electroseismogram  with 
the  scaling  used. 

In  Figure  3  the  converted  electric  and  magnetic  absolute  field  amplitudes  versus  antenna 
offsets  are  shown  at  20,  50  and  95  m  from  the  interface.  The  field  amplitude  at  each  antenna 
position  is  determined  by  calculating  a  root  mean  square  amplitude  inside  an  estimated 
pulse  window  in  the  time  domain.  The  converted  TM  wavefield,  driven  by  the  seismic  pulse 
frequency  is  diffusive,  the  real  and  imaginary  parts  of  the  TM  velocity  are  almost  equal, 
and  therefore  the  largest  signals  are  measured  by  antennas  closest  to  the  interface.  The 
frequency  content  of  the  electromagnetic  field  is  the  same  as  the  frequency  content  of  the 
incident  seismic  wave.  Since  we  are  in  the,  1,  conducting  medium  regime  we  can 

find  a  skindepth  6  =  with  //q  =  47r  x  10“^  henry /m  the  permeability  of  the  medium 

(Kong,  1990).  If  the  propagation  distance  is  much  less  than  the  skindepth,  the  near  field  of 
the  radiating  interface,  then  the  frequency  content  in  the  mechanical  and  electromagnetic 
fields  are  the  same.  The  increase  in  amplitude  with  increasing  source-antenna  offset  and  a 
later  decay  in  amplitude  with  offset  show  similarities  with  amplitudes  that  would  be  recorded 
if  the  interface  was  replaced  by  a  seismically  induced  electric  dipole  right  under  the  source 
(the  amplitude  drop  off  results  with  antenna  offset  show  that  the  effective  source  is  not  a 
simple  dipole  but  rather  a  more  complicated  multipole). 
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To  investigate  the  TM  wavefield  conversions  in  more  detail,  snapshots  in  time  are  calcu¬ 
lated  around  the  first  mechanical  contrast.  In  Figures  4  and  5,  four  successive  time  snapshots 
are  displayed,  showing  the  mechanical  wavefield  component,  the  top  figure,  the  electro¬ 
magnetic  component,  center  figure  and  the  electromagnetic  Er  component,  bottom  figure 
at  26.6  msec,  28.9  msec,  30.5  msec  and  32.0  msec  after  the  shot  respectively.  The  wavefields 
are  determined  at  120  by  60  receiver/antenna  positions.  The  receiver/ antenna  spacings  are 
1  m  in  the  horizontal  direction  and  3  m  in  the  vertical  direction.  At  t  =  26.6  msec  the 
Ricker  wavelet  (Hosken,  1988)  front  has  reached  the  interface  at  -90  m,  see  left  column  in 
Figure  4.  The  t  =  26.6  msec  snapshot  displaying  the  magnetic  field  component  of  the  TM 
mode  shows  an  amplitude  radiation  pattern  of  an  equivalent  magnetic  current  loop,  with  a 
position  centered  under  the  seismic  source  at  distance  60  m,  with  field  lines  pointing  at  one 
side  into  the  paper  and  at  the  other  side  out  of  the  paper.  The  magnetic  field  diffuses  quickly 
away  from  the  interface.  The  t  =  26.6  msec  snapshot  displaying  the  electric  field  component 
of  the  TM  mode  shows  an  amplitude  radiation  pattern  of  an  electric  dipole,  positioned  right 
under  the  seismic  source  at  a  distance  60  m.  The  largest  electric  fields  are  associated  with 
the  field  inside  the  P  wave  pulse.  At  later  times  when  a  larger  part  of  the  seismic  pulse  has 
traversed  the  interface  the  magnetic  current  loop  diameter  increases  and  the  T M  wavespeed 
differences  above  and  below  the  interface  become  visible.  The  current  system  imbalancess 
across  the  interface  change  in  direction  in  accordance  with  the  pulse  polarity,  the  magnetic 
and  electric  field  polarities  fiip  accordingly.  The  last  two  snapshots  show  complex  patterns  of 
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mechanical  field  convergences  into  electromagnetic  fields.  More  than  one  magnetic  current 
loops  appear  with  alternating  field  line  directions  and  the  inner  loop  diameter  increases. 
Also  more  electric  dipoles  with  opposite  “dipole  type  moments”  appear  there  where  the 
wavefront  passes.  The  combined  Er  and  amplitude  pattern  away  from  the  interface  has 
a  predominant  effective  electrical  dipole  character  away  from  the  interface.  In  both  cases 
the  amplitudes  decay  fast  with  traveled  distance. 

7.2.  The  Electroseismic  Conversion  at  Electrical  Contrasts 

The  explosion  source  is  again  positioned  100  m  above  an  electrical  contrast  in  the  upper 
halfspace.  Fifteen  receivers  are  positioned  symmetrically  in  a  straight  horizontal  line  at 
both  sides  of  the  source  95  m  above  the  interface  to  record  the  reflection  from  the  contrast 
and  95  m  below  the  interface  to  record  the  transmissions  through  the  contrast.  The  receiver 
spacing  is  10  m.  The  medium  input  parameters  describing  the  two  halfspaces  are  given  in 
Table  I.  The  calculated  medium  velocities  and  bulk  conductivities  based  on  the  medium 
parameters  listed  in  Table  I  are  shown  in  Table  III.  In  Figure  7  the  mechanical  displacement 
seismograms  {uz  and  components)  and  the  electroseismograms  (TM  mode  components) 
are  shown.  All  four  plots  are  seismogram  scaled.  At  the  electrical  contrast  the  P  wave 
reflects  only  as  converted  TM  waves.  Since  there  is  no  mechanical  contrast,  the  P  wave 
doesn’t  reflect  or  convert  into  shear  waves.  The  mechanical  Ur  and  Uz  seismograms  show 
the  transmitted  wave  fields  through  the  electrical  contrast  at  the  geophones  95  m  below  the 
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iiiterfacB  which  are  essentially  the  direct  P  wave  field  from  the  explosion  source  in  the  upper 
halfspace.  The  TM  wavefield  components  are  recorded  at  antennas  at  approximately  the  one 
way  P  wave  traveltime  from  source  to  interface.  Since  the  direct  wave  in  the  upper  halfspace 
is  omitted,  the  wavefield  seismogram  does  not  show  induced  electric  fields  within  the 
seismic  pulse  that  travel  with  the  compressional  wave  speed. 

In  Figure  8  the  converted  electric  and  magnetic  absolute  field  amplitudes  versus  antenna 
offsets  are  shown  at  20,  50  and  95  m  from  the  interface.  An  increase  in  amplitude  with 
increasing  source-antenna  offset  up  to  a  global  maximum  is  observed.  When  the  source- 
antenna  offset  is  further  increased  the  amplitudes  decrease  monotonically.  The  amplitude 
pattern  shows  similarities  with  amplitudes  that  would  be  recorded  if  the  interface  had  been 
replaced  by  a  seismically  induced  electric  dipole  centered  right  under  the  seismic  source. 
When  Figures  3  and  8  are  compared  the  field  amplitude  curves  have  similar  shape.  The 
amplitudes  however,  are  a  factor  10  bigger  for  the  converted  TM  mode  generated  at  the 
electrical  contrast  when  compared  to  the  converted  TM  mode  amplitudes  converted  at  the 
electrical  contrast. 

In  Figures  9  and  10  fonr  successive  snapshots  in  time  are  displayed,  showing  the  mechani¬ 
cal  Uz  wavefield  component,  top  figure,  the  electromagnetic  component,  center  figure  and 
the  electromagnetic  Er  component,  bottom  figure  at  26.6  Tfisec,  28.9  msec,  30.5  tnsec  and 
32.0  msec  after  the  shot,  respectively.  The  number  of  receivers/antennas  and  their  spacings 
in  horizontal  and  vertical  directions  are  identical  to  the  mechanical  contrast  case.  The  top 
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halfspace  has  a  salinity  of  1.0  mo///  and  the  lower  halfspace  has  a  salinity  of  0.001  mol/l. 
This  results  in  a  larger  conductivity  and  lower  electromagnetic  wavefield  velocity  in  the  up¬ 
per  halfspace  than  in  the  lower  halfspace  (see  Table  III).  Therefore  the  electric  fields  within 
the  incoming  seismic  P  wave  pulse  are  small  in  amplitude  when  compared  to  the  electric 
fields  within  the  mechanical  P  wave  in  the  lower  halfspace.  This  amplitude  behavior  can  be 
clearly  observed  in  the  electromagnetic  Et  component  snapshots.  Tracing  the  transient  TM 
component  wavefield  pattern  fronts  in  Figures  9  and  10,  the  T M  wavefield  velocity  difference 
in  both  halfspaces  is  evident.  The  complex  convergence  patterns  from  mechanical  wavefield 
into  electromagnetic  TM  fields  are  similar  to  the  patterns  observed  around  the  mechanical 
contrasts  (compare  Figures  4  and  5  with  9  and  10).  The  resultant  Er  and  amplitude 
pattern  away  from  the  interface  has  again  a  predominant  effective  electrical  dipole  character 
that  is  centered  under  the  source  at  the  contrast’s  depth.  The  amplitudes  of  the  TM  mode 
transients  decay  rapidly  with  traveled  distance. 

7.3.  Vertical  ElectroSeismic  Profiling  (VESP) 

The  converted  TM  mode  amplitude  behavior  with  distance,  discussed  when  the  surface 
electroseismic  calculations  were  discussed,  suggests  that  the  VESP  geometry  setting  is  an 
important  one.  Since  with  this  technique  the  antennas  are  positioned  close  to  the  target 
of  interest,  larger  converted  electromagnetic  signals  can  be  recorded  before  they  become 
too  attenuated  with  distance.  With  the  VESP  technique  the  electroseismic  method  can  be 
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applied  to  targets,  electrical  and/or  mechanical  contrasts,  at  greater  depths. 

The  first  VESP  example  is  a  three  layer  model  with  the  first  interface  at  300  m  depth 
and  the  second  interface  at  400  m  depth.  The  first  receiver/antenna  is  positioned  at  a 
horizontal  offset  of  50  m  and  a  depth  of  155  m.  The  receiver/antenna  depth  spacing  is  10 
m.  The  medium  description  with  its  most  important  electroseismic  medium  parameters  are 
shown  in  Figure  11.  In  Figure  12  the  numerically  calculated  VESP’s  are  shown.  The  top 
plot  is  the  mechanical  displacement  response,  and  the  center  and  bottom  plots  show  the 
Er  and  TM  mode  components,  respectively.  The  top  VESP  in  Figure  12  shows  the 
P  wave  reflections  and  converted  SV  wave  reflection  at  top  and  bottom  interface  of  the 
“sand  reservoir.”  The  direct  P  wave  in  the  top  halfspace  is  omitted.  The  other  two  vertical 
electroseismic  profiles  show  the  converted  TM  mode  components  at  the  two  contrasts.  The 
converted  TM  wavefield  components  show  up  at  all  antennas  at  virtually  the  same  time. 
The  high  electrical  conductivity  in  the  center  layer  attenuates  the  electric  field  amplitudes 
completely.  The  later  arrivals  in  the  Er  VESP  are  electric  fields  within  the  seismic  pulse 
that  travel  with  the  compressional  wave  speed.  The  component  VESP  shows  the  fast 
decay  of  the  converted  electroseismic  fields  with  distance.  The  converted  Htf,  component  of 
the  TM  mode  is  much  larger  in  amplitude  than  the  induced  magnetic  fields  inside  the  SV 
wave  pulse.  The  main  reason  is  the  close  antenna  position  to  the  target  of  interest. 

The  last  numerical  example  has  a  new  electrical  contrast  added  to  the  previous  model  as 
is  shown  in  Figure  13  .The  100  m  thick  “reservoir  sand”  is  now  divided  into  two  sands  with 
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identical  mechanical  properties  saturated  with  two  fluids  of  different  salinity.  In  Figure  1 
the  calculated  VESP’s  are  shown.  The  mechanical  displacement  response  is  identical  to  the 
previous  calculation  without  the  electrical  contrast.  But  the  TM  mode  component  VESP’s 
are  modified  with  an  extra  electroseismic  conversion  at  the  electrical  contrast. 

8.  Conclusions 

A  global  matrix  method  is  used  which  solves  the  macroscopic  equations  controlling  the 
coupled  electromagnetics  and  acoustics  of  porous  media  numerically  in  layered  poro-elastic 
media  driven  by  arbitrary  seismic  point  sources.  The  coupled  equations  decouple  for  isotropic 
porous  media  into  PSVTM  and  a  SHTE  coupled  wavefield  cases.  The  induced  current 
motion  plane  determines  to  which  electromagnetic  mode  the  mechanical  waves  are  coupled. 
Seismic  motion  which  generates  relative  flow,  induces  a  ’streaming’  electrical  current  due  to 
flow  of  double-layer  ions.  The  driving  force  for  the  relative  flow  is  a  combination  of  pressure 
gradients  set  up  by  the  peaks  and  troughs  of  a  compressional  wave  and  by  grain  accelerations. 
The  relative  flow  and  therefore  current  can  be  due  to  either  compressional  and  shear  waves, 
although  the  two  cases  are  of  significant  different  nature. 

Compressional  waves  traveling  through  homogeneous  porous  medium  cause  pressure  gra¬ 
dients  which  cause  charge  to  separate.  This  induces  within  the  seismic  pulse  a  system  of 
electric  fields  that  travel  with  the  compressional  wavespeed.  Rotational  waves  traveling 
through  a  homogeneous  porous  medium  cause  grain  accelerations  and  set  up  current  sheets. 
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This  induces  within  the  seismic  pulse  a  detectable  magnetic  fields  that  travel  with  the  ro¬ 
tational  wavespeed.  Therefore  when  the  seismic  pulse  passes  an  antenna  an  electric  field  is 
recorded  inside  the  P  wave  pulse  and  a  magnetic  field  is  recorded  inside  the  S  wave  pulse. 
The  seismic  pulse  doesn’t  radiate  electromagnetic  waves  away  from  the  pulse. 

Radiating  electromagnetic  wavefields  are  converted  from  seismic  waves  however,  when 
contrasts  in  mechanical  and/or  electrical  properties  are  traversed.  The  principle  features  of 
the  converted  electromagnetic  signals  are:  (1)  contacts  all  antennas  at  approximately  the 
same  time;  (2)  arrives  at  the  antennas  at  half  of  the  seismic  travel  time  at  normal  incidence 
reflected  P  waves;  and  (3)  changes  sign  on  opposite  sides  of  the  shot.  The  frequency  content 
of  the  converted  electromagnetic  field  has  the  same  frequency  content  of  the  driving  incident 
seismic  pulse,  as  long  as  the  propagation  distances  are  much  less  than  the  electromagnetic 
skin  depth.  Root  mean  square  converted  electromagnetic  amplitudes  versus  seismic  point 
source-antenna  offset  are  calculated  for  a  mechanical  porosity  contrast  and  an  electrical, 
salinity  contrast  at  different  depths.  The  amplitude  curves  are  similar  in  shape,  first  a 
strong  increase  in  amplitude  to  a  global  maximum  is  observed  with  increasing  antenna  offset 
and  next  a  monotonic  decrease  in  amplitude  with  a  further  increase  in  antenna  offset.  The 
amplitudes  decrease  rapidly  with  traveled  distance. 

Four  snapshots  in  time  show  the  the  generation  of  a  converted  wavefield  there  where 
the  seismic  wavefront  passes  the  interface  creating  current  imbalances  across  the  interface. 
Equivalent  sources  can  be  identified  with  the  conversion  to  electromagnetic  waves.  The 
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TM  component  amplitude  versus  offset  curves  recorded  at  some  distance  from  the  interface 
show  similarities  with  the  wavefield  that  would  be  obtained  if  the  interface  was  replaced  by 
an  equivalent  electric  dipole  positioned  right  beneath  the  seismic  source  at  the  contrast’s 
depth.  The  amplitude  versus  offset  results  and  the  electroseismograms  confirm  such  an 
effective  conversion  right  beneath  the  source.  The  VESP  modeling  shows  the  rapid  decay 
of  the  converted  electroseismic  signals  with  distance.  The  antennas  close  to  the  target  of 
interest  show  larger  amplitudes  in  the  converted  signal  than  electromagnetic  signals  inside 
the  seismic  pulse.  With  increasing  distance  from  the  contrast  seismic  pulse,  totally  dominates 
in  amplitude  the  converted  fields  in  the  electroseismograms  with  increasing  distance  from 
contrast  to  antenna. 
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APPENDIX  A.  ElectroSeismic  Field- Vector  Formalism  with 

Powerflow  Normalized  Eigenvectors 

In  this  appendix  the  electroseismic  field-vector-formalism  results  are  given,  which  are  used  to 
calculate  the  electroseismograms  in  this  paper.  The  reader  is  referred  to  Pride  and  Haartsen 
(1995)  for  the  derivation  of  these  results. 

A  linear  transformation  of  the  vector  is  performed,  through  it,  a  field-vector- 

formalism  is  obtained  in  which  a  decomposition  of  into  up  and  downgoing  fields 

is  manifest.  Let  be  the  8  by  1  vector  that  is  related  to  by  the  linear 

transformation, 

q(PSVTM)  _  j^{PSVTM)^{PSVTM) ^  (^.1) 

The  is  the  eigencolumn  matrix  of  system  matrix  =  AB)  and 
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defined  as, 


jyiPSVTM)  _  lj(rri  -TM) 

6(7’+™)] .  (A-2) 

The  —  superscript  denotes  downgoing  eigenvectors  and  the  +  superscript  denotes  upgoing 
wavefield  eigenvectors.  The  wavefield  eigenvector  is  derived  to  be, 

^(PSVTM)  ^  [U:„U^,W^,T:,^,T;,^,S,Hy,E:,f  .  (A-3) 


The  up  and  downgoing  compressional  eigenvectors,  bulk  displacement  normalized,  are  given 
in  the  following  matrix: 
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with  the  compressional  velocities  determined  by, 


2  [HM  -  C^] 


=  PbM  +  pe 


1  + 


L'^Pe 


H  -  2pfC 


PbM  +  Pe 


1  + 


L'^Pe 


H  -  2pjC 


[HM  -  C^)  (pePb 

6 

1 

to 

_ 1 

1/2 

(A-6) 
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The  +  sign  denotes  the  fast  compressional  (P/)  wavefield  velocity,  Vp^,  and  the  —  sign  denotes 
the  slow  compressional  (P,)  wavefield  velocity,  Vp,. 
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The  up  and  downgoing  shear  and  T M  wavefield  eigenvectors,  bulk  displacement  normalized. 
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are  given  in  the  following  matrix: 
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with  the  shear  and  TM  wavefield  velocities  determined  by, 
2  pB  -  p)Ipb 
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(A-8) 


The  +  sign  denotes  the  shear  {SV,  SH)  wavefield  velocity  and  the  -  sign  denotes  the 
EM  wavefield  velocity.  If  a  normalization  by  vertical  power  flow  is  preferred  the  following 
normalization  factors  for  each  wave  type  need  to  be  used. 

The  SV  and  TM  vertical  power  flux 


UJ 

(S.)  =  -jO 


V*  V  J  \  pj 


(A-9) 
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The  Pf  and  P*  vertical  power  flux 


(H  +  U  +  iDC  +  nilM) 


The  SH  TE  vertical  power  flux 


{S,)  = 


L0"‘ 


(q*  +  q) 


PePe 

.  P) 


LL^^nTlrG  +  1 


(A-10) 


(A-11) 


where  'yr  and  7^,  are  defined  as, 


1 _ 


IT  =  (A-12) 

\pB  -  4l 

7L  = 

pf  -  SJ 

To  obtain  the  normalization  with  respect  to  time-averaged  Poynting  power  perpendicular 
to  the  plane  wave  front,  we  dot  Eq  (A-9)  or  (A-10)  into  the  real  part  of  the  wave  direction 
vector,  which  is  Re[pv,Q,qv\.  The  velocity  v  determines  which  eigenvector  wave  type  is 
meant,  *  denotes  the  complex  conjugate  and  q  is  the  vertical  slownness  and  p  is  the  horizontal 
slowness  of  a  wave  type.  With  e  the  electric  permittivity  and  p  the  magnetic  permeability. 
The  electrokinetic  coupling  coefficient  L{u))  used  in  Eq  (A-4),  (A-5),  (A-7),  and  (A-8)  is 
defined  in  Pride  (1994). 


PB  -  g 
Pf  -  ^ 
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TABLE  I.  The  mechanical  contrast  and  the  electrical  contrast  medium  properties. 


Property 

Top  and  bottom 

halfspaces 

Sand 

layer 

Fresh  water 

halfspace 

Brine 

halfspace 

porosities  (<?!>[%]) 

15 

30 

15 

15 

dc  permeability  (A:[m^]) 

10-12 

10-12 

10-12 

10-12 

bulk  modulus  solid  (A:s[Pa]) 

3.6  X  IQi® 

3.6  X  101° 

3.6  X  IQi® 

3.6  X  IQi® 

bulk  modulus  fluid  {kf[Pa]) 

2.2  X  10^ 

2.2  X  10® 

2.2  X  10® 

2.2  X  10® 

frame  bulk  modulus  {kfr[Pa]) 

9.9  X  10^ 

9.9  X  10® 

9.9  X  10® 

9.9  X  10® 

frame  shear  modulus  (gfr[Pc]) 

9.0  X  10® 

9.0  X  10® 

9.0  X  10® 

9.0  X  10® 

fluid  viscosity  (?7[Pos]) 

1.0  X  10-3 

1.0  X  10-3 

1.0  X  10-3 

1.0  X  10-3 

density  solid  {ps[Kg/TrP]) 

2.7  X  103 

2.7  X  103 

2.7  X  103 

2.7  X  103 

density  fluid  (pf[Kg/rrP]) 

1.0  X  103 

1.0  X  103 

1.0  X  103 

1.0  X  103 

salinity  {C[mol/l]) 

1.0  X  10-3 

1.0 

1.0  X  10-3 

1.0 

temperature  iT[K]) 

298 

298 

298 

298 

permittivity  (k/) 

80 

80 

80 

80 

permittivity  (/c*) 

4 

4 

4 

4 

tortuosity  (ctoo) 

3 

3 

3 

3 
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TABLE  II.  Calculated  complex  wavefield  velocities  and  bulk  conductivities  used  in  the  me¬ 
chanical  contrast  model.  The  diffusive  slow  wave  and  EM  wavefield  velocities  are  listed  at 
zero  and  source  center  frequency  (200  Hz).  The  real  and  imaginary  parts  of  the  fast  and 
shear  waves  are  listed. 


Properties 

Upper  and  lower 

halfspace 

Sand  reservoir 

layer 

fast  wave  velocity  [m/s] 

(3282.62,  -0.1037) 

(3158.56,  -0.7419) 

slow  wave  velocity  [m/s] 

(21.2/146.2,  -13.3/-119.1) 

(16.7/111.6,  -10.5/-100.3) 

shear  wave[m/s] 

(1769.14,  -1.39) 

(1669.25, -1.53) 

TM  wave  velocity [m/ s] 

(318890/2.04  10®,  -201223/-2.02  10®) 

(7159/45665,  -4517/  -45326) 

conductivity  (cr[5'/m]) 

0.000388 

0.77 
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TABLE  III.  Calculated  complex  wavefield  velocities  and  bulk  conductivities  used  in  the  elec¬ 
trical  contrast  model.  The  diffusive  slow  wave  and  EM  wavefield  velocities  are  listed  at  zero 
and  source  center  frequency  (200  Hz).  The  real  and  imaginary  parts  of  the  fast  and  shear 
waves  are  listed. 


Properties 

Top  fresh  water 

saturated  halfspace 

Bottom  brine 

saturated  halfspace 

fast  wave  velocity  [m/s] 

(3282.62,  -0.1037) 

(3282.62,  -0.1037) 

slow  wave  velocity  [m/s] 

(21.2/146.2,  -13.3/-119.1) 

(21.2/146.1,  -13.3/-119.1) 

shear  wave[m/s] 

(1769.14, -1.39) 

(1769.14,  -1.39) 

TM  wave  velocity[m/s] 

(318890/2.04  10®,  -201223/-2.02  10®) 

(10124/64580,  -6389/-64100) 

conductivity  (cr[5/m]) 

0.000388 

0.385 
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"Thick"  Permeable  Sand  Layer 


fc  =  200  Hz 


Figure  1:  The  thick,  relative  to  the  mechanical  wavelength,  permeable  sand  medium  config¬ 
uration. 
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antenna  position  (m)  antenna  position  (m) 


Ex  component  TM  Hy  component  TM 


geophone  position  (m)  geophone  position  (m) 


Ux  component  Uz  component 


Figure  2:  The  mechanical  displacement  component  seismograms  and  the  TM  mode  compo 
nent  electroseismograms  calculated  for  the  thick  permeable  sand  medium  configuration 
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Time  slice  t  =  26.6  msec,  wavefield  components  Uz,  Hy  and  Ex 
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Time  slice  t  =  28.9  msec,  wavefield  components  Uz,  Hy  and  Ex 
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Figure  4;  Time  slice  snapshots  calculated  at  120  by  60  geophone/receiver  positions  around 
a  mechanical  contrast  at  t  =  26.6  msec  and  t  =  28.9  msec  after  explosion.  Top  snapshot 
is  Uz  component  of  mechanical  displacement  wavefield,  center  snapshot  is  component 
of  electromagnetic  TM  mode,  bottom  snapshot  is  Er  component  of  electromagnetic  TM 
mode.  The  complementary  coloring  of  the  electromagnetic  field  snapshots  indicate  the 
polarity  reversal  at  opposite  sides  of  the  snapshot. 
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Time  slice  t  =  30.5  msec,  wavefield  components  Uz,  Hy  and  Ex 


Time  slice  t  =  32.0  msec,  wavefield  components  Uz,  Hy  and  Ex 
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Figure  5:  Time  slice  snapshots  calculated  at  120  by  60  geophone/receiver  positions  around 
a  mechanical  contrast  at  ^  =  30.5  msec  and  t  =  32.0  msec  after  explosion.  Top  snapshot 
is  Uz  component  of  mechanical  displacement  wavefield,  center  snapshot  is  component 
of  electromagnetic  TM  mode,  bottom  snapshot  is  Er  component  of  electromagnetic  TM 
mode.  The  complementary  coloring  of  the  electromagnetic  field  snapshots  indicate  the 
polarity  reversal  at  opposite  sides  of  the  snapshot. 
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Fresh  Water  /  Brine 


fc  =  200Hz 


Figure  6:  The  fresh  water-brine  medium  configuration. 


antenna  position  (m)  antenna  position  (m) 


Ex  component  TM  Hy  component  TM 


Ux  component  Uz  component 


Figure  7:  The  mechanical  displacement  component  seismograms  and  the  TM  mode  compo¬ 
nent  electroseismograms  calculated  for  the  fresh  water  -  brine  medium  configuration. 
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Fields  from  10e-3  - 1.0  mo!/l  brine  contrast 

Ex  dipoles  at  95  m  (circles),  at  50  m  (squares),  at  20  m  (triangles)  from  interface 
Antenna  Position  (m) 
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Antenna  Position  (m) 
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Figure  8:  Converted  electric  and  magnetic  root  mean  square  wavefield  amplitudes  versus 
antenna  offsets  calculated  at  20,  50  and  95  m  from  an  electrical  contrast  generated  by  an 
explosive  point  source. 


163 


Time  slice  t  =  26.6  msec,  wavefield  components  Uz,  Hy  and  Ex 


Time  slice  t  =  28.9  msec,  wavefield  components  Uz,  Hy  and  Ex 
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Figure  9:  Time  slice  snapshots  calculated  at  120  by  60  geophone/receiver  positions  around 
an  electrical  contrast  at  t  =  26.6  msec  and  t  =  28.9  msec  after  explosion.  Top  snapshot 
is  Uz  component  of  mechanical  displacement  wavefield,  center  snapshot  is  H,f,  component 
of  electromagnetic  TM  mode,  bottom  snapshot  is  Er  component  of  electromagnetic  TM 
mode.  The  complementary  coloring  of  the  electromagnetic  field  snapshots  indicate  the 
polarity  reversal  at  opposite  sides  of  the  snapshot. 
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Time  slice  t  s  30.5  msec,  wavefield  components  Uz,  Hy  and  Ex 


0  10  20  30  40  50  60  70  80  90  100  110 


Distance  (m) 


0  10  20  30  4  0  50  60  70  80  90  100  110 


Distance  (m) 


Time  slice  t  =  32.0  msec,  wavefield  components  Uz,  Hy  and  Ex 
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Figure  10:  Time  slice  snapshots  calculated  at  120  by  60  geophone/receiver  positions  around 
an  electrical  contrast  at  t  —  30.5  msec  and  t  =  32.0  msec  after  explosion.  Top  snapshot 
is  Uz  component  of  mechanical  displacement  wavefield,  center  snapshot  is  component 
of  electromagnetic  T M  mode,  bottom  snapshot  is  Er  component  of  electromagnetic  TM 
mode.  The  complementary  coloring  of  the  electromagnetic  field  snapshots  indicate  the 
polarity  reversal  at  opposite  sides  of  the  snapshot. 


Vertical  ElectroSeismic  Profiling  (VESP) 


Figure  11:  Vertical  ElectroSeismic  Profiling  medium  configuration. 
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Ux  component  VESP 


time  in  (ms) 


Ex  component  VESP 


time  in  (ms) 


Hy  component  VESP 


time  in  (ms) 


Figure  12:  The  mechanical  displacement  component  seismogram  and  the  TM  mode  compo¬ 
nent  electroseismograms  calculated  for  a  Vertical  ElectroSeismic  Profiling  medium  con¬ 
figuration. 
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Vertical  ElectroSeismic  Profiling  (VESP) 


Figure  13;  Vertical  ElectroSeismic  Profiling  medium  configuration  with  additional  electrical 
contrast. 
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Ux  component  VESP 


Figure  14:  The  mechanical  displacement  component  seismogram  and  the  TM  mode  compo¬ 
nent  electroseismograms  calculated  for  a  Vertical  ElectroSeismic  Profiling  medium  con¬ 
figuration  with  additional  electrical  contrast  inside  the  sand  reservoir. 
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